/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math4.ode;
import org.apache.commons.math4.exception.DimensionMismatchException;
import org.apache.commons.math4.exception.MaxCountExceededException;
import org.apache.commons.math4.exception.NoBracketingException;
import org.apache.commons.math4.exception.NumberIsTooSmallException;
import org.apache.commons.math4.ode.AbstractIntegrator;
import org.apache.commons.math4.ode.AbstractParameterizable;
import org.apache.commons.math4.ode.ExpandableStatefulODE;
import org.apache.commons.math4.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math4.ode.FirstOrderIntegrator;
import org.apache.commons.math4.ode.JacobianMatrices;
import org.apache.commons.math4.ode.MainStateJacobianProvider;
import org.apache.commons.math4.ode.ParameterJacobianProvider;
import org.apache.commons.math4.ode.ParameterizedODE;
import org.apache.commons.math4.ode.UnknownParameterException;
import org.apache.commons.math4.ode.JacobianMatrices.MismatchedEquations;
import org.apache.commons.math4.ode.nonstiff.DormandPrince54Integrator;
import org.apache.commons.math4.stat.descriptive.SummaryStatistics;
import org.apache.commons.math4.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
public class JacobianMatricesTest {
@Test
public void testLowAccuracyExternalDifferentiation()
throws NumberIsTooSmallException, DimensionMismatchException,
MaxCountExceededException, NoBracketingException {
// this test does not really test JacobianMatrices,
// it only shows that WITHOUT this class, attempting to recover
// the jacobians from external differentiation on simple integration
// results with low accuracy gives very poor results. In fact,
// the curves dy/dp = g(b) when b varies from 2.88 to 3.08 are
// essentially noise.
// This test is taken from Hairer, Norsett and Wanner book
// Solving Ordinary Differential Equations I (Nonstiff problems),
// the curves dy/dp = g(b) are in figure 6.5
FirstOrderIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
double hP = 1.0e-12;
SummaryStatistics residualsP0 = new SummaryStatistics();
SummaryStatistics residualsP1 = new SummaryStatistics();
for (double b = 2.88; b < 3.08; b += 0.001) {
Brusselator brusselator = new Brusselator(b);
double[] y = { 1.3, b };
integ.integrate(brusselator, 0, y, 20.0, y);
double[] yP = { 1.3, b + hP };
integ.integrate(brusselator, 0, yP, 20.0, yP);
residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
}
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 500);
Assert.assertTrue(residualsP0.getStandardDeviation() > 30);
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 700);
Assert.assertTrue(residualsP1.getStandardDeviation() > 40);
}
@Test
public void testHighAccuracyExternalDifferentiation()
throws NumberIsTooSmallException, DimensionMismatchException,
MaxCountExceededException, NoBracketingException, UnknownParameterException {
FirstOrderIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
double hP = 1.0e-12;
SummaryStatistics residualsP0 = new SummaryStatistics();
SummaryStatistics residualsP1 = new SummaryStatistics();
for (double b = 2.88; b < 3.08; b += 0.001) {
ParamBrusselator brusselator = new ParamBrusselator(b);
double[] y = { 1.3, b };
integ.integrate(brusselator, 0, y, 20.0, y);
double[] yP = { 1.3, b + hP };
brusselator.setParameter("b", b + hP);
integ.integrate(brusselator, 0, yP, 20.0, yP);
residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
}
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 0.02);
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.03);
Assert.assertTrue(residualsP0.getStandardDeviation() > 0.003);
Assert.assertTrue(residualsP0.getStandardDeviation() < 0.004);
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 0.04);
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
Assert.assertTrue(residualsP1.getStandardDeviation() > 0.007);
Assert.assertTrue(residualsP1.getStandardDeviation() < 0.008);
}
@Test
public void testWrongParameterName() {
final String name = "an-unknown-parameter";
try {
ParamBrusselator brusselator = new ParamBrusselator(2.9);
brusselator.setParameter(name, 3.0);
Assert.fail("an exception should have been thrown");
} catch (UnknownParameterException upe) {
Assert.assertTrue(upe.getMessage().contains(name));
Assert.assertEquals(name, upe.getName());
}
}
@Test
public void testInternalDifferentiation()
throws NumberIsTooSmallException, DimensionMismatchException,
MaxCountExceededException, NoBracketingException,
UnknownParameterException, MismatchedEquations {
AbstractIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
double hP = 1.0e-12;
double hY = 1.0e-12;
SummaryStatistics residualsP0 = new SummaryStatistics();
SummaryStatistics residualsP1 = new SummaryStatistics();
for (double b = 2.88; b < 3.08; b += 0.001) {
ParamBrusselator brusselator = new ParamBrusselator(b);
brusselator.setParameter(ParamBrusselator.B, b);
double[] z = { 1.3, b };
double[][] dZdZ0 = new double[2][2];
double[] dZdP = new double[2];
JacobianMatrices jacob = new JacobianMatrices(brusselator, new double[] { hY, hY }, ParamBrusselator.B);
jacob.setParameterizedODE(brusselator);
jacob.setParameterStep(ParamBrusselator.B, hP);
jacob.setInitialParameterJacobian(ParamBrusselator.B, new double[] { 0.0, 1.0 });
ExpandableStatefulODE efode = new ExpandableStatefulODE(brusselator);
efode.setTime(0);
efode.setPrimaryState(z);
jacob.registerVariationalEquations(efode);
integ.setMaxEvaluations(5000);
integ.integrate(efode, 20.0);
jacob.getCurrentMainSetJacobian(dZdZ0);
jacob.getCurrentParameterJacobian(ParamBrusselator.B, dZdP);
// Assert.assertEquals(5000, integ.getMaxEvaluations());
// Assert.assertTrue(integ.getEvaluations() > 1500);
// Assert.assertTrue(integ.getEvaluations() < 2100);
// Assert.assertEquals(4 * integ.getEvaluations(), integ.getEvaluations());
residualsP0.addValue(dZdP[0] - brusselator.dYdP0());
residualsP1.addValue(dZdP[1] - brusselator.dYdP1());
}
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.02);
Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003);
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01);
}
@Test
public void testAnalyticalDifferentiation()
throws MaxCountExceededException, DimensionMismatchException,
NumberIsTooSmallException, NoBracketingException,
UnknownParameterException, MismatchedEquations {
AbstractIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
SummaryStatistics residualsP0 = new SummaryStatistics();
SummaryStatistics residualsP1 = new SummaryStatistics();
for (double b = 2.88; b < 3.08; b += 0.001) {
Brusselator brusselator = new Brusselator(b);
double[] z = { 1.3, b };
double[][] dZdZ0 = new double[2][2];
double[] dZdP = new double[2];
JacobianMatrices jacob = new JacobianMatrices(brusselator, Brusselator.B);
jacob.addParameterJacobianProvider(brusselator);
jacob.setInitialParameterJacobian(Brusselator.B, new double[] { 0.0, 1.0 });
ExpandableStatefulODE efode = new ExpandableStatefulODE(brusselator);
efode.setTime(0);
efode.setPrimaryState(z);
jacob.registerVariationalEquations(efode);
integ.setMaxEvaluations(5000);
integ.integrate(efode, 20.0);
jacob.getCurrentMainSetJacobian(dZdZ0);
jacob.getCurrentParameterJacobian(Brusselator.B, dZdP);
// Assert.assertEquals(5000, integ.getMaxEvaluations());
// Assert.assertTrue(integ.getEvaluations() > 350);
// Assert.assertTrue(integ.getEvaluations() < 510);
residualsP0.addValue(dZdP[0] - brusselator.dYdP0());
residualsP1.addValue(dZdP[1] - brusselator.dYdP1());
}
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.014);
Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003);
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01);
}
@Test
public void testFinalResult()
throws MaxCountExceededException, DimensionMismatchException,
NumberIsTooSmallException, NoBracketingException,
UnknownParameterException, MismatchedEquations {
AbstractIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
double[] y = new double[] { 0.0, 1.0 };
Circle circle = new Circle(y, 1.0, 1.0, 0.1);
JacobianMatrices jacob = new JacobianMatrices(circle, Circle.CX, Circle.CY, Circle.OMEGA);
jacob.addParameterJacobianProvider(circle);
jacob.setInitialMainStateJacobian(circle.exactDyDy0(0));
jacob.setInitialParameterJacobian(Circle.CX, circle.exactDyDcx(0));
jacob.setInitialParameterJacobian(Circle.CY, circle.exactDyDcy(0));
jacob.setInitialParameterJacobian(Circle.OMEGA, circle.exactDyDom(0));
ExpandableStatefulODE efode = new ExpandableStatefulODE(circle);
efode.setTime(0);
efode.setPrimaryState(y);
jacob.registerVariationalEquations(efode);
integ.setMaxEvaluations(5000);
double t = 18 * FastMath.PI;
integ.integrate(efode, t);
y = efode.getPrimaryState();
for (int i = 0; i < y.length; ++i) {
Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9);
}
double[][] dydy0 = new double[2][2];
jacob.getCurrentMainSetJacobian(dydy0);
for (int i = 0; i < dydy0.length; ++i) {
for (int j = 0; j < dydy0[i].length; ++j) {
Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9);
}
}
double[] dydcx = new double[2];
jacob.getCurrentParameterJacobian(Circle.CX, dydcx);
for (int i = 0; i < dydcx.length; ++i) {
Assert.assertEquals(circle.exactDyDcx(t)[i], dydcx[i], 1.0e-7);
}
double[] dydcy = new double[2];
jacob.getCurrentParameterJacobian(Circle.CY, dydcy);
for (int i = 0; i < dydcy.length; ++i) {
Assert.assertEquals(circle.exactDyDcy(t)[i], dydcy[i], 1.0e-7);
}
double[] dydom = new double[2];
jacob.getCurrentParameterJacobian(Circle.OMEGA, dydom);
for (int i = 0; i < dydom.length; ++i) {
Assert.assertEquals(circle.exactDyDom(t)[i], dydom[i], 1.0e-7);
}
}
@Test
public void testParameterizable()
throws MaxCountExceededException, DimensionMismatchException,
NumberIsTooSmallException, NoBracketingException,
UnknownParameterException, MismatchedEquations {
AbstractIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
double[] y = new double[] { 0.0, 1.0 };
ParameterizedCircle pcircle = new ParameterizedCircle(y, 1.0, 1.0, 0.1);
double hP = 1.0e-12;
double hY = 1.0e-12;
JacobianMatrices jacob = new JacobianMatrices(pcircle, new double[] { hY, hY },
ParameterizedCircle.CX, ParameterizedCircle.CY,
ParameterizedCircle.OMEGA);
jacob.setParameterizedODE(pcircle);
jacob.setParameterStep(ParameterizedCircle.CX, hP);
jacob.setParameterStep(ParameterizedCircle.CY, hP);
jacob.setParameterStep(ParameterizedCircle.OMEGA, hP);
jacob.setInitialMainStateJacobian(pcircle.exactDyDy0(0));
jacob.setInitialParameterJacobian(ParameterizedCircle.CX, pcircle.exactDyDcx(0));
jacob.setInitialParameterJacobian(ParameterizedCircle.CY, pcircle.exactDyDcy(0));
jacob.setInitialParameterJacobian(ParameterizedCircle.OMEGA, pcircle.exactDyDom(0));
ExpandableStatefulODE efode = new ExpandableStatefulODE(pcircle);
efode.setTime(0);
efode.setPrimaryState(y);
jacob.registerVariationalEquations(efode);
integ.setMaxEvaluations(50000);
double t = 18 * FastMath.PI;
integ.integrate(efode, t);
y = efode.getPrimaryState();
for (int i = 0; i < y.length; ++i) {
Assert.assertEquals(pcircle.exactY(t)[i], y[i], 1.0e-9);
}
double[][] dydy0 = new double[2][2];
jacob.getCurrentMainSetJacobian(dydy0);
for (int i = 0; i < dydy0.length; ++i) {
for (int j = 0; j < dydy0[i].length; ++j) {
Assert.assertEquals(pcircle.exactDyDy0(t)[i][j], dydy0[i][j], 5.0e-4);
}
}
double[] dydp0 = new double[2];
jacob.getCurrentParameterJacobian(ParameterizedCircle.CX, dydp0);
for (int i = 0; i < dydp0.length; ++i) {
Assert.assertEquals(pcircle.exactDyDcx(t)[i], dydp0[i], 5.0e-4);
}
double[] dydp1 = new double[2];
jacob.getCurrentParameterJacobian(ParameterizedCircle.OMEGA, dydp1);
for (int i = 0; i < dydp1.length; ++i) {
Assert.assertEquals(pcircle.exactDyDom(t)[i], dydp1[i], 1.0e-2);
}
}
private static class Brusselator extends AbstractParameterizable
implements MainStateJacobianProvider, ParameterJacobianProvider {
public static final String B = "b";
private double b;
public Brusselator(double b) {
super(B);
this.b = b;
}
@Override
public int getDimension() {
return 2;
}
@Override
public void computeDerivatives(double t, double[] y, double[] yDot) {
double prod = y[0] * y[0] * y[1];
yDot[0] = 1 + prod - (b + 1) * y[0];
yDot[1] = b * y[0] - prod;
}
@Override
public void computeMainStateJacobian(double t, double[] y, double[] yDot,
double[][] dFdY) {
double p = 2 * y[0] * y[1];
double y02 = y[0] * y[0];
dFdY[0][0] = p - (1 + b);
dFdY[0][1] = y02;
dFdY[1][0] = b - p;
dFdY[1][1] = -y02;
}
@Override
public void computeParameterJacobian(double t, double[] y, double[] yDot,
String paramName, double[] dFdP) {
if (isSupported(paramName)) {
dFdP[0] = -y[0];
dFdP[1] = y[0];
} else {
dFdP[0] = 0;
dFdP[1] = 0;
}
}
public double dYdP0() {
return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
}
public double dYdP1() {
return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
}
}
private static class ParamBrusselator extends AbstractParameterizable
implements FirstOrderDifferentialEquations, ParameterizedODE {
public static final String B = "b";
private double b;
public ParamBrusselator(double b) {
super(B);
this.b = b;
}
@Override
public int getDimension() {
return 2;
}
/** {@inheritDoc} */
@Override
public double getParameter(final String name)
throws UnknownParameterException {
complainIfNotSupported(name);
return b;
}
/** {@inheritDoc} */
@Override
public void setParameter(final String name, final double value)
throws UnknownParameterException {
complainIfNotSupported(name);
b = value;
}
@Override
public void computeDerivatives(double t, double[] y, double[] yDot) {
double prod = y[0] * y[0] * y[1];
yDot[0] = 1 + prod - (b + 1) * y[0];
yDot[1] = b * y[0] - prod;
}
public double dYdP0() {
return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
}
public double dYdP1() {
return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
}
}
/** ODE representing a point moving on a circle with provided center and angular rate. */
private static class Circle extends AbstractParameterizable
implements MainStateJacobianProvider, ParameterJacobianProvider {
public static final String CX = "cx";
public static final String CY = "cy";
public static final String OMEGA = "omega";
private final double[] y0;
private double cx;
private double cy;
private double omega;
public Circle(double[] y0, double cx, double cy, double omega) {
super(CX,CY,OMEGA);
this.y0 = y0.clone();
this.cx = cx;
this.cy = cy;
this.omega = omega;
}
@Override
public int getDimension() {
return 2;
}
@Override
public void computeDerivatives(double t, double[] y, double[] yDot) {
yDot[0] = omega * (cy - y[1]);
yDot[1] = omega * (y[0] - cx);
}
@Override
public void computeMainStateJacobian(double t, double[] y,
double[] yDot, double[][] dFdY) {
dFdY[0][0] = 0;
dFdY[0][1] = -omega;
dFdY[1][0] = omega;
dFdY[1][1] = 0;
}
@Override
public void computeParameterJacobian(double t, double[] y, double[] yDot,
String paramName, double[] dFdP)
throws UnknownParameterException {
complainIfNotSupported(paramName);
if (paramName.equals(CX)) {
dFdP[0] = 0;
dFdP[1] = -omega;
} else if (paramName.equals(CY)) {
dFdP[0] = omega;
dFdP[1] = 0;
} else {
dFdP[0] = cy - y[1];
dFdP[1] = y[0] - cx;
}
}
public double[] exactY(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
double dx0 = y0[0] - cx;
double dy0 = y0[1] - cy;
return new double[] {
cx + cos * dx0 - sin * dy0,
cy + sin * dx0 + cos * dy0
};
}
public double[][] exactDyDy0(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
return new double[][] {
{ cos, -sin },
{ sin, cos }
};
}
public double[] exactDyDcx(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
return new double[] {1 - cos, -sin};
}
public double[] exactDyDcy(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
return new double[] {sin, 1 - cos};
}
public double[] exactDyDom(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
double dx0 = y0[0] - cx;
double dy0 = y0[1] - cy;
return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) };
}
}
/** ODE representing a point moving on a circle with provided center and angular rate. */
private static class ParameterizedCircle extends AbstractParameterizable
implements FirstOrderDifferentialEquations, ParameterizedODE {
public static final String CX = "cx";
public static final String CY = "cy";
public static final String OMEGA = "omega";
private final double[] y0;
private double cx;
private double cy;
private double omega;
public ParameterizedCircle(double[] y0, double cx, double cy, double omega) {
super(CX,CY,OMEGA);
this.y0 = y0.clone();
this.cx = cx;
this.cy = cy;
this.omega = omega;
}
@Override
public int getDimension() {
return 2;
}
@Override
public void computeDerivatives(double t, double[] y, double[] yDot) {
yDot[0] = omega * (cy - y[1]);
yDot[1] = omega * (y[0] - cx);
}
@Override
public double getParameter(final String name)
throws UnknownParameterException {
if (name.equals(CX)) {
return cx;
} else if (name.equals(CY)) {
return cy;
} else if (name.equals(OMEGA)) {
return omega;
} else {
throw new UnknownParameterException(name);
}
}
@Override
public void setParameter(final String name, final double value)
throws UnknownParameterException {
if (name.equals(CX)) {
cx = value;
} else if (name.equals(CY)) {
cy = value;
} else if (name.equals(OMEGA)) {
omega = value;
} else {
throw new UnknownParameterException(name);
}
}
public double[] exactY(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
double dx0 = y0[0] - cx;
double dy0 = y0[1] - cy;
return new double[] {
cx + cos * dx0 - sin * dy0,
cy + sin * dx0 + cos * dy0
};
}
public double[][] exactDyDy0(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
return new double[][] {
{ cos, -sin },
{ sin, cos }
};
}
public double[] exactDyDcx(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
return new double[] {1 - cos, -sin};
}
public double[] exactDyDcy(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
return new double[] {sin, 1 - cos};
}
public double[] exactDyDom(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
double dx0 = y0[0] - cx;
double dy0 = y0[1] - cy;
return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) };
}
}
}