/**
* 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 java.util.Arrays;
import org.apache.commons.lang.Validate;
/**
*
*/
public class PDEGrid1D {
private final int _nSpaceNodes;
private final double[] _tNodes;
private final double[] _xNodes;
private final double[] _dt;
private final double[] _dx;
private final double[][] _x1st;
private final double[][] _x1stFwd;
private final double[][] _x1stBkd;
private final double[][] _x2nd;
/**
* Create a uniform grid with numTimeNodes between 0 and tMax, and numSpaceNodes between xMin and xMax
* @param numTimeNodes The number of time nodes. Note, the number of time steps is numTimeNodes - 1, therefore need numTimeNodes >= 2
* @param numSpaceNodes The number of space nodes. Note, this includes the boundaries, so the number of internal nodes is numSpaceNodes - 2, therefore need numSpaceNodes >=3
* @param tMax maximum time
* @param xMin minimum x
* @param xMax maximum x
*/
public PDEGrid1D(final int numTimeNodes, final int numSpaceNodes, final double tMax, final double xMin, final double xMax) {
Validate.isTrue(numTimeNodes > 1, "need at least 2 time nodes");
Validate.isTrue(numSpaceNodes > 2, "need at least 3 space nodes");
Validate.isTrue(tMax > 0, "need tMax > 0");
Validate.isTrue(xMax > xMin, "need xMax > xMin");
_nSpaceNodes = numSpaceNodes;
_tNodes = new double[numTimeNodes];
_xNodes = new double[numSpaceNodes];
_dt = new double[numTimeNodes - 1];
_dx = new double[numSpaceNodes - 1];
final double dt = tMax / (numTimeNodes - 1);
_tNodes[numTimeNodes - 1] = tMax;
for (int i = 0; i < numTimeNodes - 1; i++) {
_tNodes[i] = i * dt;
_dt[i] = dt;
}
final double dx = (xMax - xMin) / (numSpaceNodes - 1);
_xNodes[numSpaceNodes - 1] = xMax;
for (int i = 0; i < numSpaceNodes - 1; i++) {
_xNodes[i] = xMin + i * dx;
_dx[i] = dx;
}
_x1st = new double[numSpaceNodes - 2][3];
_x2nd = new double[numSpaceNodes - 2][3];
for (int i = 0; i < numSpaceNodes - 2; i++) {
_x1st[i][0] = -1. / 2. / dx;
_x1st[i][1] = 0.0;
_x1st[i][2] = 1. / 2. / dx;
_x2nd[i][0] = 1. / dx / dx;
_x2nd[i][1] = -2. / dx / dx;
_x2nd[i][2] = 1. / dx / dx;
}
_x1stFwd = new double[numSpaceNodes - 1][2];
_x1stBkd = new double[numSpaceNodes - 1][2];
for (int i = 0; i < numSpaceNodes - 1; i++) {
_x1stFwd[i][0] = -1 / dx;
_x1stFwd[i][1] = 1 / dx;
_x1stBkd[i][0] = -1 / dx;
_x1stBkd[i][1] = 1 / dx;
}
}
public PDEGrid1D(final MeshingFunction timeMesh, final MeshingFunction spaceMesh) {
this(timeMesh.getPoints(), spaceMesh.getPoints());
}
public PDEGrid1D(final double[] timeGrid, final double[] spaceGrid) {
final int tNodes = timeGrid.length;
final int xNodes = spaceGrid.length;
Validate.isTrue(tNodes > 1, "need at least 2 time nodes");
Validate.isTrue(xNodes > 2, "need at least 3 space nodes");
_nSpaceNodes = xNodes;
_tNodes = timeGrid;
_xNodes = spaceGrid;
_dt = new double[tNodes - 1];
for (int n = 0; n < tNodes - 1; n++) {
_dt[n] = timeGrid[n + 1] - timeGrid[n];
Validate.isTrue(_dt[n] > 0, "time steps must be increasing");
}
_dx = new double[xNodes - 1];
for (int i = 0; i < xNodes - 1; i++) {
_dx[i] = spaceGrid[i + 1] - spaceGrid[i];
Validate.isTrue(_dx[i] > 0, "space steps must be increasing");
}
_x1st = new double[xNodes - 2][3];
_x2nd = new double[xNodes - 2][3];
for (int i = 0; i < xNodes - 2; i++) {
_x1st[i][0] = -_dx[i + 1] / _dx[i] / (_dx[i] + _dx[i + 1]);
_x1st[i][1] = (_dx[i + 1] - _dx[i]) / _dx[i] / _dx[i + 1];
_x1st[i][2] = _dx[i] / _dx[i + 1] / (_dx[i] + _dx[i + 1]);
_x2nd[i][0] = 2 / _dx[i] / (_dx[i] + _dx[i + 1]);
_x2nd[i][1] = -2 / _dx[i] / _dx[i + 1];
_x2nd[i][2] = 2 / _dx[i + 1] / (_dx[i] + _dx[i + 1]);
}
_x1stFwd = new double[xNodes - 1][2];
_x1stBkd = new double[xNodes - 1][2];
for (int i = 0; i < xNodes - 1; i++) {
_x1stFwd[i][0] = -1 / _dx[i];
_x1stFwd[i][1] = 1 / _dx[i];
_x1stBkd[i][0] = -1 / _dx[i];
_x1stBkd[i][1] = 1 / _dx[i];
}
}
public int getNumTimeNodes() {
return _tNodes.length;
}
public int getNumSpaceNodes() {
return _xNodes.length;
}
public double getTimeNode(final int n) {
return _tNodes[n];
}
public double getSpaceNode(final int i) {
return _xNodes[i];
}
public double getTimeStep(final int n) {
return _dt[n];
}
public double getSpaceStep(final int i) {
return _dx[i];
}
public double[] getTimeNodes() {
return _tNodes;
}
public double[] getSpaceNodes() {
return _xNodes;
}
public int getLowerBoundIndexForTime(final double time) {
return getLowerBoundIndex(_tNodes, time);
}
public int getLowerBoundIndexForSpace(final double space) {
return getLowerBoundIndex(_xNodes, space);
}
public double[] getFirstDerivativeCoefficients(final int i) {
Validate.isTrue(i > 0 && i < _nSpaceNodes - 1, "Can't take central difference at first or last node. Use Forward or backwards");
return _x1st[i - 1];
}
public double[] getFirstDerivativeForwardCoefficients(final int i) {
Validate.isTrue(i < _nSpaceNodes - 1, "Can't take forward difference at last node. Use central or backwards");
return _x1stFwd[i];
}
public double[] getFirstDerivativeBackwardCoefficients(final int i) {
Validate.isTrue(i > 0, "Can't take backwards difference at first node. Use central or forwards");
return _x1stBkd[i - 1];
}
public double[] getSecondDerivativeCoefficients(final int i) {
if (i == 0) {
return _x2nd[0]; // TODO check this is still the best 3-point when the grid is non-uniform
} else if (i == _nSpaceNodes - 1) {
return _x2nd[_nSpaceNodes - 3];
}
return _x2nd[i - 1];
}
public PDEGrid1D withDoubleTimeSteps() {
final double[] timeGrid = new double[_tNodes.length * 2 - 1];
for (int i = 0; i < _tNodes.length - 1; i++) {
timeGrid[2 * i] = _tNodes[i];
timeGrid[2 * i + 1] = (_tNodes[i] + _tNodes[i + 1]) / 2.0;
}
timeGrid[2 * (_tNodes.length - 1)] = _tNodes[_tNodes.length - 1];
return new PDEGrid1D(timeGrid, _xNodes);
}
private int getLowerBoundIndex(final double[] array, final double t) {
final int n = array.length;
if (t < array[0]) {
return 0;
}
if (t > array[n - 1]) {
return n - 1;
}
int index = Arrays.binarySearch(array, t);
if (index >= 0) {
// Fast break out if it's an exact match.
return index;
}
if (index < 0) {
index = -(index + 1);
index--;
}
return index;
}
@Override
public int hashCode() {
final int prime = 31;
int result = 1;
result = prime * result + Arrays.hashCode(_dt);
result = prime * result + Arrays.hashCode(_dx);
result = prime * result + _nSpaceNodes;
result = prime * result + Arrays.hashCode(_tNodes);
result = prime * result + Arrays.hashCode(_x1st);
result = prime * result + Arrays.hashCode(_x1stBkd);
result = prime * result + Arrays.hashCode(_x1stFwd);
result = prime * result + Arrays.hashCode(_x2nd);
result = prime * result + Arrays.hashCode(_xNodes);
return result;
}
@Override
public boolean equals(final Object obj) {
if (this == obj) {
return true;
}
if (obj == null) {
return false;
}
if (getClass() != obj.getClass()) {
return false;
}
final PDEGrid1D other = (PDEGrid1D) obj;
if (!Arrays.equals(_dt, other._dt)) {
return false;
}
if (!Arrays.equals(_dx, other._dx)) {
return false;
}
if (_nSpaceNodes != other._nSpaceNodes) {
return false;
}
if (!Arrays.equals(_tNodes, other._tNodes)) {
return false;
}
if (!Arrays.equals(_x1st, other._x1st)) {
return false;
}
if (!Arrays.equals(_x1stBkd, other._x1stBkd)) {
return false;
}
if (!Arrays.equals(_x1stFwd, other._x1stFwd)) {
return false;
}
if (!Arrays.equals(_x2nd, other._x2nd)) {
return false;
}
if (!Arrays.equals(_xNodes, other._xNodes)) {
return false;
}
return true;
}
}