//TODO shunt this setup into its own class
ConvectionDiffusionPDE1DStandardCoefficients pde = pdeProvider.getForwardLocalVol(lvsm);
Function1D<Double, Double> initialCond = initialConProvider.getForwardCallPut(true);
double xL = 0.8;
double xH = 1.2;
BoundaryCondition lower = new NeumannBoundaryCondition(-1.0, xL, true);
BoundaryCondition upper = new NeumannBoundaryCondition(0.0, xH, false);
final MeshingFunction spaceMeshF = new HyperbolicMeshing(xL, xH, 1.0, 200, 0.001);
final MeshingFunction timeMeshF = new ExponentialMeshing(0, t, 50, 4.0);
final MeshingFunction timeMeshB = new DoubleExponentialMeshing(0, t, t / 2, 50, 2.0, -4.0);
final PDEGrid1D grid = new PDEGrid1D(timeMeshF, spaceMeshF);
PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dbF = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pde, initialCond, lower, upper, grid);
PDETerminalResults1D res = (PDETerminalResults1D) solver.solve(dbF);
final double minK = Math.exp(-6 * rootT);
final double maxK = Math.exp(6 * rootT);
Map<Double, Double> vols = PDEUtilityTools.priceToImpliedVol(fwdCurve, t, res, minK, maxK, true);
DoubleQuadraticInterpolator1D interpolator = Interpolator1DFactory.DOUBLE_QUADRATIC_INSTANCE;
Interpolator1DDataBundle idb = interpolator.getDataBundle(vols);
//set up for solving backwards PDE
ConvectionDiffusionPDE1DStandardCoefficients pdeB = pdeProvider.getBackwardsLocalVol(t, lvsm);
double sL = xL * spot;
double sH = xH * spot;
final MeshingFunction spaceMeshB = new HyperbolicMeshing(sL, sH, spot, 200, 0.001);
final PDEGrid1D gridB = new PDEGrid1D(timeMeshB, spaceMeshB);
int index = SurfaceArrayUtils.getLowerBoundIndex(gridB.getSpaceNodes(), spot);
double s1 = gridB.getSpaceNode(index);
double s2 = gridB.getSpaceNode(index + 1);
final double w = (s2 - spot) / (s2 - s1);
//solve a separate backwards PDE for each strike
for (int i = 0; i < 10; i++) {
double z = -5 + i;
double k = spot * Math.exp(0.4 * rootT * z);
double x = k / fwd;
double vol = ivs.getVolatility(t, k);
double volFPDE = interpolator.interpolate(idb, x);
boolean isCall = (k >= fwd);
BoundaryCondition lowerB = new NeumannBoundaryCondition(isCall ? 0 : -1, sL, true);
BoundaryCondition upperB = new NeumannBoundaryCondition(isCall ? 1 : 0, sH, false);
Function1D<Double, Double> bkdIC = initialConProvider.getEuropeanPayoff(k, isCall);
PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dbB = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pdeB, bkdIC, lowerB, upperB, gridB);
PDEResults1D resB = solver.solve(dbB);
double price1 = resB.getFunctionValue(index);