001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.ode; 019 020 /** 021 * This class is used in the junit tests for the ODE integrators. 022 023 * <p>This specific problem is the following differential equation : 024 * <pre> 025 * y1'' = -y1/r^3 y1 (0) = 1-e y1' (0) = 0 026 * y2'' = -y2/r^3 y2 (0) = 0 y2' (0) =sqrt((1+e)/(1-e)) 027 * r = sqrt (y1^2 + y2^2), e = 0.9 028 * </pre> 029 * This is a two-body problem in the plane which can be solved by 030 * Kepler's equation 031 * <pre> 032 * y1 (t) = ... 033 * </pre> 034 * </p> 035 036 */ 037 public class TestProblem3 038 extends TestProblemAbstract { 039 040 /** Serializable version identifier. */ 041 private static final long serialVersionUID = 8567328542728919999L; 042 043 /** Eccentricity */ 044 double e; 045 046 /** theoretical state */ 047 private double[] y; 048 049 /** 050 * Simple constructor. 051 * @param e eccentricity 052 */ 053 public TestProblem3(double e) { 054 super(); 055 this.e = e; 056 double[] y0 = { 1 - e, 0, 0, Math.sqrt((1+e)/(1-e)) }; 057 setInitialConditions(0.0, y0); 058 setFinalConditions(20.0); 059 double[] errorScale = { 1.0, 1.0, 1.0, 1.0 }; 060 setErrorScale(errorScale); 061 y = new double[y0.length]; 062 } 063 064 /** 065 * Simple constructor. 066 */ 067 public TestProblem3() { 068 this(0.1); 069 } 070 071 /** 072 * Copy constructor. 073 * @param problem problem to copy 074 */ 075 public TestProblem3(TestProblem3 problem) { 076 super(problem); 077 e = problem.e; 078 y = problem.y.clone(); 079 } 080 081 /** {@inheritDoc} */ 082 public TestProblem3 copy() { 083 return new TestProblem3(this); 084 } 085 086 @Override 087 public void doComputeDerivatives(double t, double[] y, double[] yDot) { 088 089 // current radius 090 double r2 = y[0] * y[0] + y[1] * y[1]; 091 double invR3 = 1 / (r2 * Math.sqrt(r2)); 092 093 // compute the derivatives 094 yDot[0] = y[2]; 095 yDot[1] = y[3]; 096 yDot[2] = -invR3 * y[0]; 097 yDot[3] = -invR3 * y[1]; 098 099 } 100 101 @Override 102 public double[] computeTheoreticalState(double t) { 103 104 // solve Kepler's equation 105 double E = t; 106 double d = 0; 107 double corr = 999.0; 108 for (int i = 0; (i < 50) && (Math.abs(corr) > 1.0e-12); ++i) { 109 double f2 = e * Math.sin(E); 110 double f0 = d - f2; 111 double f1 = 1 - e * Math.cos(E); 112 double f12 = f1 + f1; 113 corr = f0 * f12 / (f1 * f12 - f0 * f2); 114 d -= corr; 115 E = t + d; 116 }; 117 118 double cosE = Math.cos(E); 119 double sinE = Math.sin(E); 120 121 y[0] = cosE - e; 122 y[1] = Math.sqrt(1 - e * e) * sinE; 123 y[2] = -sinE / (1 - e * cosE); 124 y[3] = Math.sqrt(1 - e * e) * cosE / (1 - e * cosE); 125 126 return y; 127 } 128 129 }