source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/integration/RungeKutta.java

Last change on this file was 127, checked in by Wouter Pasman, 6 years ago

#41 ROLL BACK of rev.126 . So this version is equal to rev. 125

File size: 28.4 KB
Line 
1/*
2* Class RungeKutta
3* requires interfaces DerivFunction and DerivnFunction
4*
5* Contains the methods for the Runge-Kutta procedures for solving
6* single or solving sets of ordinary differential equations (ODEs)
7* [draws heavily on the approach adopted in Numerical Recipes
8* (C language version)http://www.nr.com]
9*
10* A single ODE is supplied by means of an interface,
11* DerivFunction
12* A set of ODEs is supplied by means of an interface,
13* DerivnFunction
14*
15* WRITTEN BY: Dr Michael Thomas Flanagan
16*
17* DATE: February 2002
18* UPDATES: 22 June 2003, April 2004,
19* 15 September 2006 (to incorporate improvements suggested by Klaus Benary [Klaus.Benary@gede.de])
20* 11 April 2007, 25 April 2007, 4 July 2008, 26-31 January 2010
21*
22* DOCUMENTATION:
23* See Michael Thomas Flanagan's Java library on-line web page:
24* http://www.ee.ucl.ac.uk/~mflanaga/java/RungeKutta.html
25* http://www.ee.ucl.ac.uk/~mflanaga/java/
26*
27* Copyright (c) 2002 - 2010
28*
29* PERMISSION TO COPY:
30* Permission to use, copy and modify this software and its documentation for
31* NON-COMMERCIAL purposes is granted, without fee, provided that an acknowledgement
32* to the author, Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies.
33*
34* Dr Michael Thomas Flanagan makes no representations about the suitability
35* or fitness of the software for any or for a particular purpose.
36* Michael Thomas Flanagan shall not be liable for any damages suffered
37* as a result of using, modifying or distributing this software or its derivatives.
38*
39***************************************************************************************/
40
41package agents.anac.y2015.agentBuyogV2.flanagan.integration;
42
43import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
44
45
46// Class for Runge-Kutta solution of ordinary differential equations
47public class RungeKutta{
48
49 private double x0 = Double.NaN; // initial value of x
50 private double xn = Double.NaN;; // final value of x
51 private double y0 = Double.NaN; // initial value of y; single ODE
52 private double[] yy0 = null; // initial values of y; multiple ODEs
53 private int nODE = 0; // number of ODEs
54 private double step = Double.NaN; // step size
55
56 private double relTol = 1.0e-5; // tolerance multiplicative factor in adaptive step methods
57 private double absTol = 1.0e-3; // tolerance additive factor to ensure non-zeto tolerance in adaptive step methods
58 private int maxIter = -1; // maximum iterations allowed in adaptive step methods
59 private int nIter = 0; // number of iterations taken
60 private static int nStepsMultiplier = 1000; // multiplied by number of steps to give internally calculated
61 // maximum allowed iterations in adaptive step methods
62
63 private static double safetyFactor = 0.9; // safety factor for Runge Kutta Fehlberg and CashKarp tolerance checks as
64 // error estimations are uncertain
65 private static double incrementFactor = -0.2; // factor used in calculating a step size increment in Fehlberg and CashKarp procedures
66 private static double decrementFactor = -0.25; // factor used in calculating a step size decrement in Fehlberg and CashKarp procedures
67
68 public RungeKutta(){
69 }
70
71 public void setInitialValueOfX(double x0){
72 this.x0 = x0;
73 }
74
75 public void setFinalValueOfX(double xn){
76 this.xn = xn;
77 }
78
79 public void setInitialValueOfY(double y0){
80 this.setInitialValuesOfY(y0);
81 }
82
83 public void setInitialValueOfY(double[] yy0){
84 this.setInitialValuesOfY(yy0);
85 }
86
87 public void setInitialValuesOfY(double y0){
88 this.y0 = y0;
89 this.yy0 = new double[1];
90 this.yy0[0] = y0;
91 this.nODE = 1;
92 }
93
94 public void setInitialValuesOfY(double[] yy0){
95 this.yy0 = yy0;
96 this.nODE = yy0.length;
97 if(this.nODE==1)this.y0 = yy0[0];
98 }
99
100 public void setStepSize(double step){
101 this.step = step;
102 }
103
104 public void setToleranceScalingFactor(double relTol){
105 this.relTol = relTol;
106 }
107
108 public void setToleranceAdditionFactor(double absTol){
109 this.absTol = absTol;
110 }
111
112 public void setMaximumIterations(int maxIter){
113 this.maxIter = maxIter;
114 }
115
116 public int getNumberOfIterations(){
117 return this.nIter;
118 }
119
120 public static void resetNstepsMultiplier(int multiplier){
121 RungeKutta.nStepsMultiplier = multiplier;
122 }
123
124 // Fourth order Runge-Kutta for a single ordinary differential equation
125 // Non-static method
126 public double fourthOrder(DerivFunction g){
127 if(Double.isNaN(this.x0))throw new IllegalArgumentException("No initial x value has been entered");
128 if(Double.isNaN(this.xn))throw new IllegalArgumentException("No final x value has been entered");
129 if(Double.isNaN(this.y0))throw new IllegalArgumentException("No initial y value has been entered");
130 if(Double.isNaN(this.step))throw new IllegalArgumentException("No step size has been entered");
131
132 double k1 = 0.0D, k2 = 0.0D, k3 = 0.0D, k4 = 0.0D;
133 double x = 0.0D, y = this.y0;
134
135 // Calculate nsteps
136 double ns = (this.xn - this.x0)/this.step;
137 ns = Math.rint(ns);
138 int nsteps = (int) ns; // number of steps
139 this.nIter = nsteps;
140 double stepUsed = (this.xn - this.x0)/ns;
141
142 for(int i=0; i<nsteps; i++){
143 x = this.x0 + i*stepUsed;
144
145 k1 = stepUsed*g.deriv(x, y);
146 k2 = stepUsed*g.deriv(x + stepUsed/2, y + k1/2);
147 k3 = stepUsed*g.deriv(x + stepUsed/2, y + k2/2);
148 k4 = stepUsed*g.deriv(x + stepUsed, y + k3);
149
150 y += k1/6 + k2/3 + k3/3 + k4/6;
151 }
152 return y;
153 }
154
155
156 // Fourth order Runge-Kutta for a single ordinary differential equation (ODE)
157 // Static method
158 public static double fourthOrder(DerivFunction g, double x0, double y0, double xn, double h){
159
160 RungeKutta rk = new RungeKutta();
161 rk.setInitialValueOfX(x0);
162 rk.setFinalValueOfX(xn);
163 rk.setInitialValueOfY(y0);
164 rk.setStepSize(h);
165
166 return rk.fourthOrder(g);
167 }
168
169 // Fourth order Runge-Kutta for n (nODE) ordinary differential equations (ODE)
170 // Non-static method
171 public double[] fourthOrder(DerivnFunction g){
172 if(Double.isNaN(this.x0))throw new IllegalArgumentException("No initial x value has been entered");
173 if(Double.isNaN(this.xn))throw new IllegalArgumentException("No final x value has been entered");
174 if(this.yy0==null)throw new IllegalArgumentException("No initial y values have been entered");
175 if(Double.isNaN(this.step))throw new IllegalArgumentException("No step size has been entered");
176
177 double[] k1 =new double[this.nODE];
178 double[] k2 =new double[this.nODE];
179 double[] k3 =new double[this.nODE];
180 double[] k4 =new double[this.nODE];
181 double[] y =new double[this.nODE];
182 double[] yd =new double[this.nODE];
183 double[] dydx =new double[this.nODE];
184 double x = 0.0D;
185
186 // Calculate nsteps
187 double ns = (this.xn - this.x0)/this.step;
188 ns = Math.rint(ns);
189 int nsteps = (int) ns;
190 this.nIter = nsteps;
191 double stepUsed = (this.xn - this.x0)/ns;
192
193 // initialise
194 for(int i=0; i<this.nODE; i++)y[i] = this.yy0[i];
195
196 // iteration over allowed steps
197 for(int j=0; j<nsteps; j++){
198 x = this.x0 + j*stepUsed;
199 dydx = g.derivn(x, y);
200 for(int i=0; i<this.nODE; i++)k1[i] = stepUsed*dydx[i];
201
202 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + k1[i]/2;
203 dydx = g.derivn(x + stepUsed/2, yd);
204 for(int i=0; i<this.nODE; i++)k2[i] = stepUsed*dydx[i];
205
206 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + k2[i]/2;
207 dydx = g.derivn(x + stepUsed/2, yd);
208 for(int i=0; i<this.nODE; i++)k3[i] = stepUsed*dydx[i];
209
210 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + k3[i];
211 dydx = g.derivn(x + stepUsed, yd);
212 for(int i=0; i<this.nODE; i++)k4[i] = stepUsed*dydx[i];
213
214 for(int i=0; i<this.nODE; i++)y[i] += k1[i]/6 + k2[i]/3 + k3[i]/3 + k4[i]/6;
215
216 }
217 return y;
218 }
219
220 // Fourth order Runge-Kutta for n ordinary differential equations (ODE)
221 // Static method
222 public static double[] fourthOrder(DerivnFunction g, double x0, double[] y0, double xn, double h){
223
224 RungeKutta rk = new RungeKutta();
225 rk.setInitialValueOfX(x0);
226 rk.setFinalValueOfX(xn);
227 rk.setInitialValuesOfY(y0);
228 rk.setStepSize(h);
229
230 return rk.fourthOrder(g);
231 }
232
233
234 // Runge-Kutta-Cash-Karp for a single ordinary differential equation (ODE)
235 // Non-static method
236 public double cashKarp(DerivFunction g){
237 if(Double.isNaN(this.x0))throw new IllegalArgumentException("No initial x value has been entered");
238 if(Double.isNaN(this.xn))throw new IllegalArgumentException("No final x value has been entered");
239 if(Double.isNaN(this.y0))throw new IllegalArgumentException("No initial y value has been entered");
240 if(Double.isNaN(this.step))throw new IllegalArgumentException("No step size has been entered");
241
242 double k1 = 0.0D, k2 = 0.0D, k3 = 0.0D, k4 = 0.0D, k5 = 0.0D, k6 = 0.0D;
243 double y = this.y0, y5 = 0.0D, y6 = 0.0D, yd = 0.0D, dydx = 0.0D;
244 double x = this.x0, err = 0.0D, delta = 0.0D, tol = 0.0D;
245 int counter = 0;
246
247 if(this.maxIter==-1){
248 this.maxIter = (int)(RungeKutta.nStepsMultiplier*(this.xn - this.x0)/this.step);
249 }
250 double stepUsed = this.step;
251
252
253 while(x<this.xn){
254 counter++;
255 if(counter>this.maxIter)throw new ArithmeticException("Maximum number of iterations exceeded");
256 dydx = g.deriv(x, y);
257 k1 = stepUsed*dydx;
258
259 yd = y + k1/5.0;
260 dydx = g.deriv(x + stepUsed/5.0, yd);
261 k2 = stepUsed*dydx;
262
263 yd = y + (3.0*k1 + 9.0*k2)/40.0;
264 dydx = g.deriv(x + 3.0*stepUsed/10.0, yd);
265 k3 = stepUsed*dydx;
266
267 yd = y + (3.0*k1 - 9.0*k2 + 12.0*k3)/10.0;
268 dydx = g.deriv(x + 3.0*stepUsed/5.0, yd);
269 k4 = stepUsed*dydx;
270
271 yd = y -11.0*k1/54.0 + 5.0*k2/2.0 - 70.0*k3/27.0 + 35.0*k4/27.0;
272 dydx = g.deriv(x + stepUsed, yd);
273 k5 = stepUsed*dydx;
274
275 yd = y + 1631.0*k1/55296.0 + 175.0*k2/512.0 + 575.0*k3/13824.0 + 44275.0*k4/110592.0 + 253.0*k5/4096.0;
276 dydx = g.deriv(x + 7.0*stepUsed/8.0, yd);
277 k6 = stepUsed*dydx;
278
279 y5 = y + 2825.0*k1/27648.0 + 18575.0*k3/48384.0 + 13525.0*k4/55296.0 + 277.0*k5/14336.0 + k6/4.0;
280 y6 = y + 37*k1/378.0 + 250.0*k3/621.0 + 125.0*k4/594.0 + 512.0*k6/1771.0;
281 err = Math.abs(y6 - y5);
282 tol= err/(Math.abs(y5)*this.relTol + this.absTol);
283 if(tol<=1.0){
284 x += stepUsed;
285 delta = RungeKutta.safetyFactor*Math.pow(tol, RungeKutta.incrementFactor);
286 if(delta>4.0){
287 stepUsed*= 4.0;
288 }else if(delta >1.0){
289 stepUsed*=delta;
290 }
291 if(x+stepUsed > this.xn)stepUsed = this.xn - x;
292 y = y5;
293 }
294 else{
295 delta = RungeKutta.safetyFactor*Math.pow(tol,RungeKutta.decrementFactor);
296 if(delta < 0.1){
297 stepUsed *= 0.1;
298 }
299 else{
300 stepUsed *= delta;
301 }
302 }
303 }
304 this.nIter = counter;
305 return y;
306 }
307
308 // Runge-Kutta-Cash-Karp for a single ordinary differential equation (ODE)
309 // Static method
310 public static double cashKarp(DerivFunction g, double x0, double y0, double xn, double h, double absTol, double relTol, int maxIter){
311
312 RungeKutta rk = new RungeKutta();
313 rk.setInitialValueOfX(x0);
314 rk.setFinalValueOfX(xn);
315 rk.setInitialValueOfY(y0);
316 rk.setStepSize(h);
317 rk.setToleranceScalingFactor(relTol);
318 rk.setToleranceAdditionFactor(absTol);
319 rk.setMaximumIterations(maxIter);
320
321 return rk.cashKarp(g);
322 }
323
324 // Runge-Kutta-Cash-Karp for a single ordinary differential equation (ODE)
325 // Static method
326 // maximum iteration default option
327 public static double cashKarp(DerivFunction g, double x0, double y0, double xn, double h, double absTol, double relTol){
328
329 int maxIter = (int)(RungeKutta.nStepsMultiplier*(xn - x0)/h);
330
331 RungeKutta rk = new RungeKutta();
332 rk.setInitialValueOfX(x0);
333 rk.setFinalValueOfX(xn);
334 rk.setInitialValueOfY(y0);
335 rk.setStepSize(h);
336 rk.setToleranceScalingFactor(relTol);
337 rk.setToleranceAdditionFactor(absTol);
338 rk.setMaximumIterations(maxIter);
339
340 return rk.cashKarp(g);
341 }
342
343 // Runge-Kutta-Cash-Karp for n ordinary differential equations (ODEs)
344 // Non-static method
345 public double[] cashKarp(DerivnFunction g){
346 if(Double.isNaN(this.x0))throw new IllegalArgumentException("No initial x value has been entered");
347 if(Double.isNaN(this.xn))throw new IllegalArgumentException("No final x value has been entered");
348 if(this.yy0==null)throw new IllegalArgumentException("No initial y values have been entered");
349 if(Double.isNaN(this.step))throw new IllegalArgumentException("No step size has been entered");
350
351 double[] k1 =new double[this.nODE];
352 double[] k2 =new double[this.nODE];
353 double[] k3 =new double[this.nODE];
354 double[] k4 =new double[this.nODE];
355 double[] k5 =new double[this.nODE];
356 double[] k6 =new double[this.nODE];
357 double[] y =new double[this.nODE];
358 double[] y6 =new double[this.nODE];
359 double[] y5 =new double[this.nODE];
360 double[] yd =new double[this.nODE];
361 double[] dydx =new double[this.nODE];
362
363 double err = 0.0D, maxerr = 0.0D, delta = 0.0D, tol = 1.0D;
364 int counter = 0;
365
366 // initialise
367 for(int i=0; i<this.nODE; i++)y[i] = this.yy0[i];
368 double x = this.x0;
369 if(this.maxIter==-1){
370 this.maxIter = (int)(RungeKutta.nStepsMultiplier*(this.xn - this.x0)/this.step);
371 }
372 double stepUsed = this.step;
373
374 while(x<this.xn){
375 counter++;
376 if(counter>this.maxIter)throw new ArithmeticException("Maximum number of iterations exceeded");
377
378 dydx = g.derivn(x, y);
379 for(int i=0; i<this.nODE; i++)k1[i] = stepUsed*dydx[i];
380
381 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + k1[i]/5.0;
382 dydx = g.derivn(x + stepUsed/5.0, yd);
383 for(int i=0; i<this.nODE; i++)k2[i] = stepUsed*dydx[i];
384
385 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + (3.0*k1[i] + 9.0*k2[i])/40.0;
386 dydx = g.derivn(x + 3.0*stepUsed/10.0, yd);
387 for(int i=0; i<this.nODE; i++)k3[i] = stepUsed*dydx[i];
388
389 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + (3.0*k1[i] - 9.0*k2[i] + 12.0*k3[i])/10.0;
390 dydx = g.derivn(x + 3.0*stepUsed/5.0, yd);
391 for(int i=0; i<this.nODE; i++)k4[i] = stepUsed*dydx[i];
392
393 for(int i=0; i<this.nODE; i++)yd[i] = y[i] -11.0*k1[i]/54.0 + 5.0*k2[i]/2.0 - 70.0*k3[i]/27.0 + 35.0*k4[i]/27.0;
394 dydx = g.derivn(x + stepUsed, yd);
395 for(int i=0; i<this.nODE; i++)k5[i] = stepUsed*dydx[i];
396
397 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + 1631.0*k1[i]/55296.0 + 175.0*k2[i]/512.0 + 575.0*k3[i]/13824.0 + 44275.0*k4[i]/110592.0 + 253.0*k5[i]/4096.0;
398 dydx = g.derivn(x + 7.0*stepUsed/8.0, yd);
399 for(int i=0; i<this.nODE; i++)k6[i] = stepUsed*dydx[i];
400
401 maxerr=0.0D;
402 for(int i=0; i<this.nODE; i++){
403 y5[i] = y[i] + 2825.0*k1[i]/27648.0 + 18575.0*k3[i]/48384.0 + 13525.0*k4[i]/55296.0 + 277.0*k5[i]/14336.0 + k6[i]/4.0;
404 y6[i] = y[i] + 37*k1[i]/378.0 + 250.0*k3[i]/621.0 + 125.0*k4[i]/594.0 + 512.0*k6[i]/1771.0;
405 err = Math.abs(y6[i] - y5[i]);
406 tol= Math.abs(y5[i])*this.relTol + this.absTol;
407 maxerr = Math.max(maxerr,err/tol);
408 }
409 if(maxerr<=1.0D){
410 x += stepUsed;
411 delta = RungeKutta.safetyFactor*Math.pow(maxerr, RungeKutta.incrementFactor);
412 if(delta>4.0){
413 stepUsed*= 4.0;
414 }
415 else if (delta > 1.0){
416 stepUsed*=delta;
417 }
418 if(x+stepUsed > this.xn)stepUsed = this.xn-x;
419 y = Conv.copy(y5);
420 }
421 else{
422 delta = RungeKutta.safetyFactor*Math.pow(maxerr,RungeKutta.decrementFactor);
423 if(delta < 0.1D){
424 stepUsed *= 0.1;
425 }
426 else{
427 stepUsed *= delta;
428 }
429 }
430 }
431 this.nIter = counter;
432 return y;
433 }
434
435 // Runge-Kutta-Cash-Karp for n ordinary differential equations (ODEs)
436 // Static method
437 public static double[] cashKarp(DerivnFunction g, double x0, double[] y0, double xn, double h, double absTol, double relTol, int maxIter){
438
439 RungeKutta rk = new RungeKutta();
440 rk.setInitialValueOfX(x0);
441 rk.setFinalValueOfX(xn);
442 rk.setInitialValuesOfY(y0);
443 rk.setStepSize(h);
444 rk.setToleranceScalingFactor(relTol);
445 rk.setToleranceAdditionFactor(absTol);
446 rk.setMaximumIterations(maxIter);
447
448 return rk.cashKarp(g);
449 }
450
451 public static double[] cashKarp(DerivnFunction g, double x0, double[] y0, double xn, double h, double absTol, double relTol){
452
453 double nsteps = (xn - x0)/h;
454 int maxIter = (int) nsteps*RungeKutta.nStepsMultiplier;
455
456 return cashKarp(g, x0, y0, xn, h, absTol, relTol, maxIter);
457 }
458
459 // Runge-Kutta-Fehlberg for a single ordinary differential equation (ODE)
460 // Non-static method
461 public double fehlberg(DerivFunction g){
462 if(Double.isNaN(this.x0))throw new IllegalArgumentException("No initial x value has been entered");
463 if(Double.isNaN(this.xn))throw new IllegalArgumentException("No final x value has been entered");
464 if(Double.isNaN(this.y0))throw new IllegalArgumentException("No initial y value has been entered");
465 if(Double.isNaN(this.step))throw new IllegalArgumentException("No step size has been entered");
466
467 double k1 = 0.0D, k2 = 0.0D, k3 = 0.0D, k4 = 0.0D, k5 = 0.0D, k6 = 0.0D;
468 double x = this.x0, y = this.y0, y5 = 0.0D, y6 = 0.0D, err = 0.0D, delta = 0.0D, tol = 0.0D;
469 int counter = 0;
470
471 if(this.maxIter==-1){
472 this.maxIter = (int)(RungeKutta.nStepsMultiplier*(this.xn - this.x0)/this.step);
473 }
474 double stepUsed = this.step;
475
476 while(x<this.xn){
477 counter++;
478 if(counter>this.maxIter)throw new ArithmeticException("Maximum number of iterations exceeded");
479 k1 = stepUsed*g.deriv(x, y);
480 k2 = stepUsed*g.deriv(x + stepUsed/4.0, y + k1/4.0);
481 k3 = stepUsed*g.deriv(x + 3.0*stepUsed/8.0, y + (3.0*k1 + 9.0*k2)/32.0);
482 k4 = stepUsed*g.deriv(x + 12.0*stepUsed/13.0, y + (1932.0*k1 - 7200.0*k2 + 7296.0*k3)/2197.0);
483 k5 = stepUsed*g.deriv(x + stepUsed, y + 439.0*k1/216.0 - 8.0*k2 + 3680.0*k3/513.0 - 845*k4/4104.0);
484 k6 = stepUsed*g.deriv(x + 0.5*stepUsed, y - 8.0*k1/27.0 + 2.0*k2 - 3544.0*k3/2565.0 + 1859.0*k4/4104.0 - 11.0*k5/40.0);
485
486 y5 = y + 25.0*k1/216.0 + 1408.0*k3/2565.0 + 2197.0*k4/4104.0 - k5/5.0;
487 y6 = y + 16.0*k1/135.0 + 6656.0*k3/12825.0 + 28561.0*k4/56430.0 - 9.0*k5/50.0 + 2.0*k6/55.0;
488 err = Math.abs(y6 - y5);
489 tol= err/(Math.abs(y5)*this.relTol + this.absTol);
490 if(tol<=1.0){
491 x += stepUsed;
492 delta = RungeKutta.safetyFactor*Math.pow(tol, RungeKutta.incrementFactor);
493 if(delta>4.0){
494 stepUsed*= 4.0;
495 }else if(delta <1.0){
496 stepUsed*=delta;
497 }
498 if(x+stepUsed > this.xn)stepUsed = this.xn-x;
499 y = y5;
500 }
501 else{
502 delta = RungeKutta.safetyFactor*Math.pow(tol,RungeKutta.decrementFactor);
503 if(delta < 0.1){
504 stepUsed *= 0.1;
505 }
506 else{
507 stepUsed *= delta;
508 }
509 }
510 }
511 this.nIter = counter;
512 return y;
513 }
514
515 // Runge-Kutta-Fehlberg for a single ordinary differential equation (ODE)
516 // Static method
517 public static double fehlberg(DerivFunction g, double x0, double y0, double xn, double h, double absTol, double relTol, int maxIter){
518 RungeKutta rk = new RungeKutta();
519 rk.setInitialValueOfX(x0);
520 rk.setFinalValueOfX(xn);
521 rk.setInitialValueOfY(y0);
522 rk.setStepSize(h);
523 rk.setToleranceScalingFactor(relTol);
524 rk.setToleranceAdditionFactor(absTol);
525 rk.setMaximumIterations(maxIter);
526
527 return rk.fehlberg(g);
528 }
529
530 // Runge-Kutta-Fehlberg for a single ordinary differential equation (ODE)
531 // Maximum iteration default option
532 // Static method
533 public static double fehlberg(DerivFunction g, double x0, double y0, double xn, double h, double absTol, double relTol){
534
535 double nsteps = (xn - x0)/h;
536 int maxIter = (int) nsteps*RungeKutta.nStepsMultiplier;
537
538 return fehlberg(g, x0, y0, xn, h, absTol, relTol, maxIter);
539 }
540
541 // Runge-Kutta-Fehlberg for n ordinary differential equations (ODEs)
542 // Non-static method
543 public double[] fehlberg(DerivnFunction g){
544 if(Double.isNaN(this.x0))throw new IllegalArgumentException("No initial x value has been entered");
545 if(Double.isNaN(this.xn))throw new IllegalArgumentException("No final x value has been entered");
546 if(this.yy0==null)throw new IllegalArgumentException("No initial y values have been entered");
547 if(Double.isNaN(this.step))throw new IllegalArgumentException("No step size has been entered");
548
549 double[] k1 =new double[this.nODE];
550 double[] k2 =new double[this.nODE];
551 double[] k3 =new double[this.nODE];
552 double[] k4 =new double[this.nODE];
553 double[] k5 =new double[this.nODE];
554 double[] k6 =new double[this.nODE];
555 double[] y =new double[this.nODE];
556 double[] y6 =new double[this.nODE];
557 double[] y5 =new double[this.nODE];
558 double[] yd =new double[this.nODE];
559 double[] dydx =new double[this.nODE];
560
561 double err = 0.0D, maxerr = 0.0D, delta = 0.0D, tol = 1.0D;
562 int counter = 0;
563
564 // initialise
565 for(int i=0; i<this.nODE; i++)y[i] = this.yy0[i];
566 double x = this.x0;
567 if(this.maxIter==-1){
568 this.maxIter = (int)(RungeKutta.nStepsMultiplier*(this.xn - this.x0)/this.step);
569 }
570 double stepUsed = this.step;
571
572 while(x<this.xn){
573 counter++;
574 if(counter>maxIter)throw new ArithmeticException("Maximum number of iterations exceeded");
575 dydx = g.derivn(x, y);
576 for(int i=0; i<this.nODE; i++)k1[i] = stepUsed*dydx[i];
577
578 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + k1[i]/4.0;
579 dydx = g.derivn(x + stepUsed/4.0, yd);
580 for(int i=0; i<this.nODE; i++)k2[i] = stepUsed*dydx[i];
581
582 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + (3.0*k1[i] + 9.0*k2[i])/32.0;
583 dydx = g.derivn(x + 3.0*stepUsed/8.0, yd);
584 for(int i=0; i<this.nODE; i++)k3[i] = stepUsed*dydx[i];
585
586 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + (1932.0*k1[i] - 7200.0*k2[i] + 7296.0*k3[i])/2197.0;
587 dydx = g.derivn(x + 12.0*stepUsed/13.0, yd);
588 for(int i=0; i<this.nODE; i++)k4[i] = stepUsed*dydx[i];
589
590 for(int i=0; i<this.nODE; i++)yd[i] = y[i] + 439.0*k1[i]/216.0 - 8.0*k2[i] + 3680.0*k3[i]/513.0 - 845*k4[i]/4104.0;
591 dydx = g.derivn(x + stepUsed, yd);
592 for(int i=0; i<this.nODE; i++)k5[i] = stepUsed*dydx[i];
593
594 for(int i=0; i<this.nODE; i++)yd[i] = y[i] - 8.0*k1[i]/27.0 + 2.0*k2[i] - 3544.0*k3[i]/2565.0 + 1859.0*k4[i]/4104.0 - 11.0*k5[i]/40.0;
595 dydx = g.derivn(x + 0.5*stepUsed, yd);
596 for(int i=0; i<this.nODE; i++)k6[i] = stepUsed*dydx[i];
597
598 maxerr=0.0D;
599 for(int i=0; i<this.nODE; i++){
600 y5[i] = y[i] + 25.0*k1[i]/216.0 + 1408.0*k3[i]/2565.0 + 2197.0*k4[i]/4104.0 - k5[i]/5.0;
601 y6[i] = y[i] + 16.0*k1[i]/135.0 + 6656.0*k3[i]/12825.0 + 28561.0*k4[i]/56430.0 - 9.0*k5[i]/50.0 + 2.0*k6[i]/55.0;
602 err = Math.abs(y6[i] - y5[i]);
603 tol= y5[i]*this.relTol + this.absTol;
604 maxerr = Math.max(maxerr,err/tol);
605 }
606
607 if(maxerr<=1.0D){
608 x += stepUsed;
609 delta = RungeKutta.safetyFactor*Math.pow(maxerr, RungeKutta.incrementFactor);
610 if(delta>4.0){
611 stepUsed*= 4.0;
612 }
613 else if(delta > 1.0){
614 stepUsed*=delta;
615 }
616 if(x+stepUsed > this.xn)stepUsed = this.xn-x;
617 y = Conv.copy(y5);
618 }
619 else{
620 delta = RungeKutta.safetyFactor*Math.pow(maxerr,RungeKutta.decrementFactor);
621 if(delta < 0.1){
622 stepUsed *= 0.1;
623 }
624 else{
625 stepUsed *= delta;
626 }
627 }
628 }
629 this.nIter = counter;
630 return y;
631 }
632
633 // Runge-Kutta-Fehlberg for n (this.nODE) ordinary differential equations (ODEs)
634 // Static method
635 public static double[] fehlberg(DerivnFunction g, double x0, double[] y0, double xn, double h, double absTol, double relTol, int maxIter){
636
637 RungeKutta rk = new RungeKutta();
638 rk.setInitialValueOfX(x0);
639 rk.setFinalValueOfX(xn);
640 rk.setInitialValuesOfY(y0);
641 rk.setStepSize(h);
642 rk.setToleranceScalingFactor(relTol);
643 rk.setToleranceAdditionFactor(absTol);
644 rk.setMaximumIterations(maxIter);
645
646 return rk.fehlberg(g);
647 }
648
649 // Runge-Kutta-Fehlberg for n (this.nODE) ordinary differential equations (ODEs)
650 // Static method
651 // maximum iteration default option
652 public static double[] fehlberg(DerivnFunction g, double x0, double[] y0, double xn, double h, double absTol, double relTol){
653
654 double nsteps = (xn - x0)/h;
655 int maxIter = (int) nsteps*RungeKutta.nStepsMultiplier;
656
657 return fehlberg(g, x0, y0, xn, h, absTol, relTol, maxIter);
658 }
659}
660
661
662
663
Note: See TracBrowser for help on using the repository browser.