com.opengamma.analytics.financial.model.finitedifference.CraigSneydFiniteDifference2D.java Source code

Java tutorial

Introduction

Here is the source code for com.opengamma.analytics.financial.model.finitedifference.CraigSneydFiniteDifference2D.java

Source

/**
 * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
 *
 * Please see distribution for license.
 */
package com.opengamma.analytics.financial.model.finitedifference;

import org.apache.commons.lang.Validate;

import com.opengamma.analytics.math.cube.Cube;

/**
 * Craig-Sneyd splitting
 * <b>Note</b> this is for testing purposes and is not recommended for actual use
 *
 */
@SuppressWarnings("deprecation")
public class CraigSneydFiniteDifference2D implements ConvectionDiffusionPDESolver2D {

    // private static final Decomposition<?> DCOMP = new LUDecompositionCommons();
    // Theta = 0 - explicit
    private static final double THETA = 0.5;

    @Override
    public double[][] solve(final ConvectionDiffusion2DPDEDataBundle pdeData, final int tSteps, final int xSteps,
            final int ySteps, final double tMax, final BoundaryCondition2D xLowerBoundary,
            final BoundaryCondition2D xUpperBoundary, final BoundaryCondition2D yLowerBoundary,
            final BoundaryCondition2D yUpperBoundary) {
        return solve(pdeData, tSteps, xSteps, ySteps, tMax, xLowerBoundary, xUpperBoundary, yLowerBoundary,
                yUpperBoundary, null);
    }

    @Override
    public double[][] solve(final ConvectionDiffusion2DPDEDataBundle pdeData, final int tSteps, final int xSteps,
            final int ySteps, final double tMax, final BoundaryCondition2D xLowerBoundary,
            final BoundaryCondition2D xUpperBoundary, final BoundaryCondition2D yLowerBoundary,
            final BoundaryCondition2D yUpperBoundary, final Cube<Double, Double, Double, Double> freeBoundary) {

        final double dt = tMax / (tSteps);
        final double dx = (xUpperBoundary.getLevel() - xLowerBoundary.getLevel()) / (xSteps);
        final double dy = (yUpperBoundary.getLevel() - yLowerBoundary.getLevel()) / (ySteps);
        final double dtdx2 = dt / dx / dx;
        final double dtdx = dt / dx;
        final double dtdy2 = dt / dy / dy;
        final double dtdy = dt / dy;
        final double dtdxdy = dt / dx / dy;

        final double[][] v = new double[xSteps + 1][ySteps + 1];

        final double[][] vt = new double[xSteps + 1][ySteps + 1];
        final double[] x = new double[xSteps + 1];
        final double[] y = new double[ySteps + 1];

        final double[] q = new double[xSteps + 1];
        final double[] r = new double[ySteps + 1];
        final double[][] mx = new double[xSteps + 1][xSteps + 1];
        final double[][] my = new double[ySteps + 1][ySteps + 1];
        initializeMatrices(pdeData, xSteps, ySteps, xLowerBoundary, yLowerBoundary, dx, dy, v, x, y);

        double t = 0.0;
        double a, b, c, d, e, f;

        for (int n = 0; n < tSteps; n++) {

            // stag 1 full Explicit
            for (int i = 1; i < xSteps; i++) {
                for (int j = 1; j < ySteps; j++) {
                    a = pdeData.getA(t, x[i], y[j]);
                    b = pdeData.getB(t, x[i], y[j]);
                    c = pdeData.getC(t, x[i], y[j]);
                    d = pdeData.getD(t, x[i], y[j]);
                    e = pdeData.getE(t, x[i], y[j]);
                    f = pdeData.getF(t, x[i], y[j]);

                    vt[i][j] = (1 - dt * (1 - 0.5 * THETA) * c) * v[i][j];
                    vt[i][j] -= dtdx2 * a * (1 - THETA) * (v[i + 1][j] + v[i - 1][j] - 2 * v[i][j]);
                    vt[i][j] -= 0.5 * dtdx * b * (1 - THETA) * (v[i + 1][j] - v[i - 1][j]);
                    vt[i][j] -= dtdy2 * d * (v[i][j + 1] + v[i][j - 1] - 2 * v[i][j]);
                    // upwind
                    // if (f > 0) {
                    // vt[i][j] -= dtdy * f * (v[i][j] - v[i][j - 1]);
                    // } else if (f < 0) {
                    // vt[i][j] -= dtdy * f * (v[i][j + 1] - v[i][j]);
                    // }
                    vt[i][j] -= 0.5 * dtdy * f * (v[i][j + 1] - v[i][j - 1]);
                    vt[i][j] -= 0.25 * dtdxdy * e
                            * (v[i + 1][j + 1] + v[i - 1][j - 1] - v[i + 1][j - 1] - v[i - 1][j + 1]);
                }

                // really not sure what to do with boundary conditions in these intermediate steps
                vt[i][0] = v[i][0];
                vt[i][ySteps] = v[i][ySteps];
            }

            // for (int i = 0; i <= xSteps; i++) {
            // double[] temp = yLowerBoundary.getRightMatrixCondition(pdeData, t, x[i]);
            // double sum = 0;
            // for (int k = 0; k < temp.length; k++) {
            // sum += temp[k] * v[i][k];
            // }
            // sum += yLowerBoundary.getConstant(pdeData, t, x[i], dy);
            //
            // temp = yLowerBoundary.getLeftMatrixCondition(pdeData, t, x[i]);
            // for (int k = 1; k < temp.length; k++) {
            // sum -= temp[k] * vt[i][k];
            // }
            // vt[i][0] = sum / temp[0];
            //
            // temp = yUpperBoundary.getRightMatrixCondition(pdeData, t, x[i]);
            // sum = 0;
            // for (int k = 0; k < temp.length; k++) {
            // sum += temp[k] * v[i][ySteps - k];
            // }
            // sum += yUpperBoundary.getConstant(pdeData, t, x[i], dy);
            //
            // temp = yUpperBoundary.getLeftMatrixCondition(pdeData, t, x[i]);
            // for (int k = 1; k < temp.length; k++) {
            // sum -= temp[k] * vt[i][ySteps - k];
            // }
            // vt[i][ySteps] = sum / temp[0];
            // }
            //
            // for (int j = 1; j < ySteps; j++) {
            // double[] temp = xLowerBoundary.getRightMatrixCondition(pdeData, t, y[j]);
            // double sum = 0;
            // for (int k = 0; k < temp.length; k++) {
            // sum += temp[k] * v[k][j];
            // }
            // sum += xLowerBoundary.getConstant(pdeData, t, y[j], dx);
            //
            // temp = xLowerBoundary.getLeftMatrixCondition(pdeData, t, y[j]);
            // for (int k = 1; k < temp.length; k++) {
            // sum -= temp[k] * vt[k][j];
            // }
            // vt[0][j] = sum / temp[0];
            //
            // temp = xUpperBoundary.getRightMatrixCondition(pdeData, t, y[j]);
            // sum = 0;
            // for (int k = 0; k < temp.length; k++) {
            // sum += temp[k] * v[xSteps - k][j];
            // }
            // sum += xUpperBoundary.getConstant(pdeData, t, y[j], dx);
            //
            // temp = xUpperBoundary.getLeftMatrixCondition(pdeData, t, y[j]);
            // for (int k = 1; k < temp.length; k++) {
            // sum -= temp[k] * vt[xSteps - k][j];
            // }
            // vt[xSteps][j] = sum / temp[0];
            // }

            // stag 2 implicit in x
            t += dt / 2;
            for (int j = 0; j <= ySteps; j++) {
                for (int i = 1; i < xSteps; i++) {
                    a = pdeData.getA(t, x[i], y[j]);
                    b = pdeData.getB(t, x[i], y[j]);
                    c = pdeData.getC(t, x[i], y[j]);

                    mx[i][i - 1] = THETA * (dtdx2 * a - 0.5 * dtdx * b);
                    mx[i][i] = 1 + THETA * (-2 * dtdx2 * a + 0.5 * dt * c);
                    mx[i][i + 1] = THETA * (dtdx2 * a + 0.5 * dtdx * b);

                    q[i] = vt[i][j];
                }

                double[] temp = xLowerBoundary.getLeftMatrixCondition(t, y[j]);
                for (int k = 0; k < temp.length; k++) {
                    mx[0][k] = temp[k];
                }
                temp = xUpperBoundary.getLeftMatrixCondition(t, y[j]);
                for (int k = 0; k < temp.length; k++) {
                    mx[xSteps][xSteps - k] = temp[k];
                }

                temp = xLowerBoundary.getRightMatrixCondition(t, y[j]);
                double sum = 0;
                for (int k = 0; k < temp.length; k++) {
                    sum += temp[k] * v[k][j];
                }
                q[0] = sum + xLowerBoundary.getConstant(t, y[j], dx);

                temp = xUpperBoundary.getRightMatrixCondition(t, y[j]);
                sum = 0;
                for (int k = 0; k < temp.length; k++) {
                    sum += temp[k] * v[xSteps - k][j];
                }
                q[xSteps] = sum + xUpperBoundary.getConstant(t, y[j], dx);

                // SOR
                final double omega = 1.5;
                final int count = sor(xSteps, vt, q, mx, j, omega);
                Validate.isTrue(count < 1000, "SOR exceeded max iterations");
            }

            // stag 3 explicit in y
            for (int i = 0; i <= xSteps; i++) {
                for (int j = 1; j < ySteps; j++) {

                    c = pdeData.getC(t, x[i], y[j]);
                    d = pdeData.getD(t, x[i], y[j]);
                    f = pdeData.getF(t, x[i], y[j]);

                    vt[i][j] += THETA * 0.5 * dt * c * v[i][j];
                    vt[i][j] += THETA * dtdy2 * d * (v[i][j + 1] + v[i][j - 1] - 2 * v[i][j]);

                    // upwind
                    // if (f > 0) {
                    // vt[i][j] += THETA * dtdy * f * (v[i][j] - v[i][j - 1]);
                    // } else if (f < 0) {
                    // vt[i][j] += THETA * dtdy * f * (v[i][j + 1] - v[i][j]);
                    // }
                    vt[i][j] += THETA * 0.5 * dtdy * f * (v[i][j + 1] - v[i][j - 1]);
                }

                // double[] temp = yLowerBoundary.getRightMatrixCondition(pdeData, t, x[i]);
                // double sum = 0;
                // for (int k = 0; k < temp.length; k++) {
                // sum += temp[k] * v[i][k];
                // }
                // sum += yLowerBoundary.getConstant(pdeData, t, x[i], dy);
                //
                // temp = yLowerBoundary.getLeftMatrixCondition(pdeData, t, x[i]);
                // for (int k = 1; k < temp.length; k++) {
                // sum -= temp[k] * vt[i][k];
                // }
                // vt[i][0] = sum / temp[0];
                //
                // temp = yUpperBoundary.getRightMatrixCondition(pdeData, t, x[i]);
                // sum = 0;
                // for (int k = 0; k < temp.length; k++) {
                // sum += temp[k] * v[i][ySteps - k];
                // }
                // sum += yUpperBoundary.getConstant(pdeData, t, x[i], dy);
                //
                // temp = yUpperBoundary.getLeftMatrixCondition(pdeData, t, x[i]);
                // for (int k = 1; k < temp.length; k++) {
                // sum -= temp[k] * vt[i][ySteps - k];
                // }
                // vt[i][ySteps] = sum / temp[0];
            }

            // The y = 0 and y = yStep boundary values are assumed the same as the previous sub-step
            // Again we could apply the y boundary conditions here

            // stag 4 implicit in y
            for (int i = 0; i <= xSteps; i++) {
                for (int j = 1; j < ySteps; j++) {

                    c = pdeData.getC(t, x[i], y[j]);
                    d = pdeData.getD(t, x[i], y[j]);
                    f = pdeData.getF(t, x[i], y[j]);

                    // upwind
                    // if (f > 0) {
                    // my[j][j - 1] = THETA * (dtdy2 * d - dtdy * f);
                    // my[j][j] = 1 + THETA * (-2 * dtdy2 * d + dtdy * f + 0.5 * dt * c);
                    // my[j][j + 1] = THETA * (dtdy2 * d);
                    // } else if (f < 0) {
                    // my[j][j - 1] = THETA * (dtdy2 * d);
                    // my[j][j] = 1 + THETA * (-2 * dtdy2 * d - dtdy * f + 0.5 * dt * c);
                    // my[j][j + 1] = THETA * (dtdy2 * d + dtdy * f);
                    // }
                    my[j][j - 1] = THETA * (dtdy2 * d - 0.5 * dtdy * f);
                    my[j][j] = 1 + THETA * (-2 * dtdy2 * d + 0.5 * dt * c);
                    my[j][j + 1] = THETA * (dtdy2 * d + 0.5 * dtdy * f);

                    r[j] = vt[i][j];
                }

                double[] temp = yLowerBoundary.getLeftMatrixCondition(t, x[i]);
                for (int k = 0; k < temp.length; k++) {
                    my[0][k] = temp[k];
                }
                temp = yUpperBoundary.getLeftMatrixCondition(t, x[i]);
                for (int k = 0; k < temp.length; k++) {
                    my[ySteps][ySteps - k] = temp[k];
                }

                temp = yLowerBoundary.getRightMatrixCondition(t, x[i]);
                double sum = 0;
                for (int k = 0; k < temp.length; k++) {
                    sum += temp[k] * v[i][k];
                }
                r[0] = sum + yLowerBoundary.getConstant(t, x[i], dy);

                temp = yUpperBoundary.getRightMatrixCondition(t, x[i]);
                sum = 0;
                for (int k = 0; k < temp.length; k++) {
                    sum += temp[k] * v[i][ySteps - k];
                }
                r[ySteps] = sum + yUpperBoundary.getConstant(t, x[i], dy);

                // SOR
                final double omega = 1.5;
                double scale = 1.0;
                double errorSqr = Double.POSITIVE_INFINITY;
                int count = 0;
                while (errorSqr / (scale + 1e-10) > 1e-18 && count < 1000) {
                    errorSqr = 0.0;
                    scale = 0.0;
                    int min, max;
                    for (int l = 0; l <= ySteps; l++) {
                        min = (l == ySteps ? 0 : Math.max(0, l - 1));
                        max = (l == 0 ? ySteps : Math.min(ySteps, l + 1));
                        sum = 0;
                        // for (int k = 0; k <= ySteps; k++) {
                        for (int k = min; k <= max; k++) {
                            sum += my[l][k] * v[i][k];
                        }
                        final double correction = omega / my[l][l] * (r[l] - sum);
                        // if (freeBoundary != null) {
                        // correction = Math.max(correction, freeBoundary.getZValue(t, x[j]) - f[j]);
                        // }
                        errorSqr += correction * correction;
                        v[i][l] += correction;
                        scale += v[i][l] * v[i][l];
                    }
                    count++;
                }
                Validate.isTrue(count < 1000, "SOR exceeded max interations");
            }

        } // time loop
        return v;

    }

    // TODO there is almost identical code lines 297-319
    private int sor(final int steps, final double[][] v, final double[] q, final double[][] mx, final int j,
            final double omega) {
        double sum;
        int min;
        int max;
        int count = 0;
        double scale = 1.0;
        double errorSqr = Double.POSITIVE_INFINITY;
        while (errorSqr / (scale + 1e-10) > 1e-18 && count < 1000) {
            errorSqr = 0.0;
            scale = 0.0;
            for (int l = 0; l <= steps; l++) {
                min = (l == steps ? 0 : Math.max(0, l - 1));
                max = (l == 0 ? steps : Math.min(steps, l + 1));
                sum = 0;
                // for (int k = 0; k <= xSteps; k++) {
                for (int k = min; k <= max; k++) { // mx is tri-diagonal so only need 3 steps here
                    sum += mx[l][k] * v[k][j];
                }
                final double correction = omega / mx[l][l] * (q[l] - sum);
                // if (freeBoundary != null) {
                // correction = Math.max(correction, freeBoundary.getZValue(t, x[j]) - f[j]);
                // }
                errorSqr += correction * correction;
                v[l][j] += correction;
                scale += v[l][j] * v[l][j];
            }
            count++;
        }
        return count;
    }

    private void initializeMatrices(final ConvectionDiffusion2DPDEDataBundle pdeData, final int xSteps,
            final int ySteps, final BoundaryCondition2D xLowerBoundary, final BoundaryCondition2D yLowerBoundary,
            final double dx, final double dy, final double[][] v, final double[] x, final double[] y) {
        double currentX = 0;
        double currentY = 0;

        for (int j = 0; j <= ySteps; j++) {
            currentY = yLowerBoundary.getLevel() + j * dy;
            y[j] = currentY;
        }
        for (int i = 0; i <= xSteps; i++) {
            currentX = xLowerBoundary.getLevel() + i * dx;
            x[i] = currentX;
            for (int j = 0; j <= ySteps; j++) {
                v[i][j] = pdeData.getInitialValue(x[i], y[j]);
            }
        }
    }

    // private double[][] solveSOR(double[][] m, double[][] v)
}