← Back to branch summary

~ubuntu-branches/ubuntu/trusty/commons-math/trusty

« back to all changes in this revision

Viewing changes to src/test/org/apache/commons/math/ode/TestProblem3.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-03-15 20:20:21 UTC
  • Revision ID: james.westby@ubuntu.com-20090315202021-zto3nmvqgcf3ami4
Tags: upstream-1.2
Import upstream version 1.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/test/org/apache/commons/math/ode/TestProblem3.java
 
1
/*
 
2
 * Licensed to the Apache Software Foundation (ASF) under one or more
 
3
 * contributor license agreements.  See the NOTICE file distributed with
 
4
 * this work for additional information regarding copyright ownership.
 
5
 * The ASF licenses this file to You under the Apache License, Version 2.0
 
6
 * (the "License"); you may not use this file except in compliance with
 
7
 * the License.  You may obtain a copy of the License at
 
8
 *
 
9
 *      http://www.apache.org/licenses/LICENSE-2.0
 
10
 *
 
11
 * Unless required by applicable law or agreed to in writing, software
 
12
 * distributed under the License is distributed on an "AS IS" BASIS,
 
13
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 
14
 * See the License for the specific language governing permissions and
 
15
 * limitations under the License.
 
16
 */
 
17
 
 
18
package org.apache.commons.math.ode;
 
19
 
 
20
/**
 
21
 * This class is used in the junit tests for the ODE integrators.
 
22
 
 
23
 * <p>This specific problem is the following differential equation :
 
24
 * <pre>
 
25
 *    y1'' = -y1/r^3  y1 (0) = 1-e  y1' (0) = 0
 
26
 *    y2'' = -y2/r^3  y2 (0) = 0    y2' (0) =sqrt((1+e)/(1-e))
 
27
 *    r = sqrt (y1^2 + y2^2), e = 0.9
 
28
 * </pre>
 
29
 * This is a two-body problem in the plane which can be solved by
 
30
 * Kepler's equation
 
31
 * <pre>
 
32
 *   y1 (t) = ...
 
33
 * </pre>
 
34
 * </p>
 
35
 
 
36
 */
 
37
class TestProblem3
 
38
  extends TestProblemAbstract {
 
39
 
 
40
  /** Eccentricity */
 
41
  double e;
 
42
 
 
43
  /** theoretical state */
 
44
  private double[] y;
 
45
 
 
46
  /**
 
47
   * Simple constructor.
 
48
   * @param e eccentricity
 
49
   */
 
50
  public TestProblem3(double e) {
 
51
    super();
 
52
    this.e = e;
 
53
    double[] y0 = { 1 - e, 0, 0, Math.sqrt((1+e)/(1-e)) };
 
54
    setInitialConditions(0.0, y0);
 
55
    setFinalConditions(20.0);
 
56
    double[] errorScale = { 1.0, 1.0, 1.0, 1.0 };
 
57
    setErrorScale(errorScale);
 
58
    y = new double[y0.length];
 
59
  }
 
60
 
 
61
  /**
 
62
   * Simple constructor.
 
63
   */
 
64
  public TestProblem3() {
 
65
    this(0.1);
 
66
  }
 
67
 
 
68
  /**
 
69
   * Copy constructor.
 
70
   * @param problem problem to copy
 
71
   */
 
72
  public TestProblem3(TestProblem3 problem) {
 
73
    super(problem);
 
74
    e = problem.e;
 
75
    y = (double[]) problem.y.clone();
 
76
  }
 
77
 
 
78
  /**
 
79
   * Clone operation.
 
80
   * @return a copy of the instance
 
81
   */
 
82
  public Object clone() {
 
83
    return new TestProblem3(this);
 
84
  }
 
85
 
 
86
  public void doComputeDerivatives(double t, double[] y, double[] yDot) {
 
87
 
 
88
    // current radius
 
89
    double r2 = y[0] * y[0] + y[1] * y[1];
 
90
    double invR3 = 1 / (r2 * Math.sqrt(r2));
 
91
 
 
92
    // compute the derivatives
 
93
    yDot[0] = y[2];
 
94
    yDot[1] = y[3];
 
95
    yDot[2] = -invR3  * y[0];
 
96
    yDot[3] = -invR3  * y[1];
 
97
 
 
98
  }
 
99
 
 
100
  public double[] computeTheoreticalState(double t) {
 
101
 
 
102
    // solve Kepler's equation
 
103
    double E = t;
 
104
    double d = 0;
 
105
    double corr = 0;
 
106
    do {
 
107
      double f2  = e * Math.sin(E);
 
108
      double f0  = d - f2;
 
109
      double f1  = 1 - e * Math.cos(E);
 
110
      double f12 = f1 + f1;
 
111
      corr  = f0 * f12 / (f1 * f12 - f0 * f2);
 
112
      d -= corr;
 
113
      E = t + d;
 
114
    } while (Math.abs(corr) > 1.0e-12);
 
115
 
 
116
    double cosE = Math.cos(E);
 
117
    double sinE = Math.sin(E);
 
118
 
 
119
    y[0] = cosE - e;
 
120
    y[1] = Math.sqrt(1 - e * e) * sinE;
 
121
    y[2] = -sinE / (1 - e * cosE);
 
122
    y[3] = Math.sqrt(1 - e * e) * cosE / (1 - e * cosE);
 
123
 
 
124
    return y;
 
125
  }
 
126
 
 
127
}