1 /*
2 * $Id: rkutta.i,v 1.1 2005-09-18 22:06:07 dhmunro Exp $
3 * 4th order Runge-Kutta integrator (rk_integrate, rkutta)
4 * Also Bulirsch-Stoer integrator (bs_integrate, bstoer)
5 * After routines in Numerical Recipes by Press et.al.
6 */
7 /* Copyright (c) 2005, The Regents of the University of California.
8 * All rights reserved.
9 * This file is part of yorick (http://yorick.sourceforge.net).
10 * Read the accompanying LICENSE file for details.
11 */
12
13 /* ------------------------------------------------------------------------ */
14
rk_integrate(derivative,y1,rk_x,epsilon,dx1)15 func rk_integrate(derivative, y1, rk_x, epsilon, dx1)
16 /* DOCUMENT y= rk_integrate(derivative, y1, x, epsilon, dx1)
17 integrates dydx= DERIVATIVE(y,x) beginning at (X(1),Y1) and
18 going to X(0) with fractional error EPSILON. The result is
19 the value of y at each value in the list X. If non-nil, DX1
20 will be used as initial guess for the first step size.
21 Otherwise, X(2)-X(1) will be the first step size guess.
22
23 The list of X values must be monotone -- strictly increasing
24 or strictly decreasing; the Runge-Kutta step sizes are selected
25 adaptively until the next X value would be passed, when the
26 step size is adjusted to complete the step exactly.
27
28 The external variable rk_maxits (default 10000) is the
29 maximum number of steps rk_integrate will take.
30
31 If a function rk_yscale(y,dydx,x,dx) exists, it is used
32 to compute an appropriate yscale to give the EPSILON error
33 criterion meaning. Otherwise, yscale is taken to be:
34 abs(y)+abs(dydx*dx)+1.e-30
35
36 Based on odeint from Numerical Recipes (Press, et.al.).
37 If the function you are trying to integrate is very
38 smooth, and your X values are fairly far apart, bs_integrate
39 may work better than rk_integrate.
40
41 SEE ALSO: rkutta, bs_integrate, rk_maxits, rk_minstep,
42 rk_maxstep, rk_ngood, rk_nbad, rkdumb, rk4
43 */
44 {
45 if (numberof(rk_x)<2) return array(y1, dimsof(rk_x));
46 local rk_y;
47 if (is_void(dx1)) dx1= rk_x(2)-rk_x(1);
48 rk_nstore= -1;
49 rkutta, derivative, y1, rk_x(1), rk_x(0), epsilon, dx1;
50 return rk_y;
51 }
52
rkutta(derivative,y0,x0,x1,epsilon,dx0)53 func rkutta(derivative, y0,x0, x1,epsilon, dx0)
54 /* DOCUMENT y1= rkutta(derivative, y0,x0, x1,epsilon, dx0)
55 integrates dydx= DERIVATIVE(y,x) beginning at (X0,Y0) and
56 going to X1 with fractional error EPSILON. The result is
57 the value of y at X1. DX0 will be used as the initial guess
58 for a step size.
59
60 If the external variable rk_nstore is >0, rk_y and rk_x
61 will contain up to rk_nstore intermediate values after the
62 call to rkutta. Consider using rk_integrate if you need
63 this feature; using rk_nstore gives you the results at
64 intermediate values which will tend to be closer where
65 the Runge-Kutta step size was smaller, while rk_integrate
66 forces you to specify precisely which x values you want.
67
68 The external variable rk_maxits (default 10000) is the
69 maximum number of steps rkutta will take. The variable
70 rk_minstep (default 0.0) is the minimum step size. The
71 variable rk_maxstep (default 1.e35) is the maximum step
72 size, which you may need if you are storing intermediate
73 values (particularly with bstoer).
74
75 If a function rk_yscale(y,dydx,x,dx) exists, it is used
76 to compute an appropriate yscale to give the EPSILON error
77 criterion meaning. Otherwise, yscale is taken to be:
78 abs(y)+abs(dydx*dx)+1.e-30
79
80 Based on odeint from Numerical Recipes (Press, et.al.).
81 If the function you are trying to integrate is very
82 smooth, bstoer will probably work better than rkutta.
83
84 SEE ALSO: rk_integrate, bstoer, rk_nstore, rk_maxits,
85 rk_minstep, rk_maxstep, rk_ngood, rk_nbad, rkdumb, rk4
86 */
87 {
88 extern rk_x, rk_y, rk_ngood, rk_nbad;
89 rk_ngood= rk_nbad= 0;
90 if (rk_nstore > 0) {
91 if (rk_nstore<2) rk_nstore= 2;
92 rk_x= array(double(x0), rk_nstore);
93 rk_y= array(0.+y0, rk_nstore);
94 s= [1, 1, 1]; /* see rk_store function */
95 } else if (rk_nstore < 0) {
96 x0= rk_x(1);
97 x1= rk_x(0);
98 if (anyof(rk_x(dif)*(x1-x0) <= 0.0))
99 error, "given rk_x must be monotonic";
100 rk_y= array(0.+y0, numberof(rk_x));
101 s= 2;
102 }
103
104 dxsign= sign(x1-x0);
105 dx= double(abs(dx0)*dxsign);
106 x= double(x0);
107 y= 0.+y0;
108 for (n=1 ; n<=rk_maxits ; ++n) {
109 dydx= derivative(y, x);
110 if (!is_void(rk_yscale)) yscale= rk_yscale(y,dydx,x,dx);
111 else yscale= abs(y)+abs(dydx*dx)+1.e-30;
112 if (abs(dx) > rk_maxstep) dx= dxsign*rk_maxstep;
113 if (dxsign*(x+dx-x1) > 0.0) dx= x1-x;
114 if (rk_nstore<0 &&
115 dxsign*(x+dx-rk_x(s)) > 0.0) dx= rk_x(s)-x;
116 local dxdid, dxnxt;
117 y= rkqc(y,dydx, x,dx, derivative, epsilon,yscale, dxdid,dxnxt);
118 x+= dxdid;
119 if (dxdid == dx) ++rk_ngood;
120 else ++rk_nbad;
121 if (rk_nstore>0) s= rk_store(y,x,s);
122 else if (rk_nstore<0 &&
123 dxsign*(x-rk_x(s))>=0.0) rk_y(..,s++)= y;
124 all_done= (dxsign*(x-x1) >= 0.0);
125 if (all_done) break;
126 if (abs(dxnxt) < abs(rk_minstep))
127 error, "required step less than rk_minstep";
128 dx= dxnxt;
129 }
130
131 if (rk_nstore>0) {
132 if (rk_x(s(3)) != x) {
133 s(2)= 1; /* always store final value */
134 s= rk_store(y,x,s);
135 }
136 rk_y= rk_y(..,1:s(3));
137 rk_x= rk_x(1:s(3));
138 }
139 if (!all_done) error, "exceeded rk_maxits iterations";
140 return y;
141 }
142
143 local rk_nstore, rk_maxits, rk_minstep, rk_maxstep, rk_ngood, rk_nbad;
144 /* DOCUMENT rk_nstore, rk_maxits, rk_minstep, rk_maxstep,
145 rk_ngood, rk_nbad
146
147 rk_nstore maximum number of y values rkutta (bstoer) will store
148 after rkutta (bstoer) call, rk_y and rk_x contain stored values
149
150 The other variables are inputs or outputs for rkutta, bstoer,
151 rk_integrate, or bs_integrate:
152
153 rk_maxits maximum number of steps (default 10000)
154 rk_minstep minimum step size (default 0.0)
155 rk_maxstep maximum step size (default 1.e35)
156 rk_ngood number of good steps taken
157 rk_nbad number of failed (but repaired) steps taken
158 */
159 rk_maxits= 10000;
160 rk_minstep= 0.0;
161 rk_maxstep= 1.0e35;
162 rk_nstore= 0;
163
rk_store(y,x,s)164 func rk_store(y,x,s)
165 {
166 /* s= [step number, step increment, store index] */
167 i= ++s(1);
168 if (! ((i-1)%s(2))) {
169 i= ++s(3);
170 if (i > rk_nstore) {
171 y2= rk_y(..,1:0:2);
172 x2= rk_x(1:0:2);
173 i= numberof(x2);
174 rk_y(..,1:i)= y2;
175 rk_x(1:i)= x2;
176 s(3)= ++i;
177 s(2)*= 2;
178 }
179 rk_y(..,i)= y;
180 rk_x(i)= x;
181 }
182 return s;
183 }
184
185 func rkdumb(derivative, y0,x0, x1,nsteps, nosave=)
186 /* DOCUMENT y_of_x= rkdumb(derivative, y0,x0, x1,nsteps)
187 integrates dydx= DERIVATIVE(y,x) beginning at (X0,Y0) and
188 going to X1 in NSTEPS 4th order Runge-Kutta steps. The
189 result is dimsof(Y0)-by-(NSTEPS+1) values of y at the points
190 span(X0, X1, NSTEPS+1).
191 If the nosave= keyword is non-zero, the returned value will
192 simply be the final y value.
193 */
194 {
195 dx= (x1-x0)/nsteps;
196 ++nsteps;
197 if (!nosave) y= array(0.+y0, nsteps);
198 for (i=2 ; i<=nsteps ; ++i) {
199 y0= rk4(y0,derivative(y0,x0), x0,dx, derivative);
200 x0+= dx;
201 if (!nosave) y(..,i)= y0;
202 }
203 return nosave? y0 : y;
204 }
205
206 func rkqc(y,dydx, x,dx, derivative, epsilon,yscale, &dxdid,&dxnxt)
207 {
208 x0= x;
209 y0= y;
210
211 for (;;) {
212 x= x0+dx;
213 if (x==x0) error, "integration step crash";
214 /* first take two half steps... */
215 dx2= 0.5*dx;
216 x2= x0+dx2;
217 y2= rk4(y0,dydx, x0,dx2, derivative);
218 y2= rk4(y2,derivative(y2,x2), x2,dx2, derivative);
219 /* ...then compare with one full step... */
220 y1= rk4(y0,dydx, x0,dx, derivative);
221 /* ...to estimate error */
222 y1= y2-y1;
223 err= max(abs(y1/yscale))/epsilon;
224 if (err <= 1.0) {
225 dxdid= dx;
226 dxnxt= (err>6.e-4)? 0.9*dx*err^-0.20 : 4.*dx;
227 break;
228 }
229 dx*= 0.9*err^-0.25;
230 }
231 return y2 + y1/15.;
232 }
233
rk4(y,dydx,x,dx,derivative)234 func rk4(y,dydx, x,dx, derivative)
235 /* DOCUMENT y_at_x_plus_dx= rk4(y,dydx, x,dx, derivative)
236 takes a single 4th order Runge-Kutta step from X to X+DX.
237 DERIVATIVE(y,x) is a function returning dydx; the input DYDX
238 is DERIVATIVE(y,x) at the input (X,Y). This fourth evaluation
239 of DERIVATIVE must be performed by the caller of rk4.
240 */
241 {
242 dx2= 0.5*dx;
243 x2= x+dx2;
244 dydxp= derivative(y+dydx*dx2, x2); /* slope at 1st trial midpoint */
245 dydxm= derivative(y+dydxp*dx2, x2); /* slope at 2nd trial midpoint */
246 dydxp+= dydxm;
247 dydxm= derivative(y+dydxm*dx, x+dx); /* slope at trial endpoint */
248 return y + (dydx+dydxp+dydxp+dydxm)*(dx/6.0);
249 }
250
251 /* ------------------------------------------------------------------------ */
252
bs_integrate(derivative,y1,rk_x,epsilon,dx1)253 func bs_integrate(derivative, y1, rk_x, epsilon, dx1)
254 /* DOCUMENT y= bs_integrate(derivative, y1, x, epsilon, dx1)
255 Bulirsch-Stoer integrator, otherwise identical to rk_integrate
256 routine. All of the options for rk_integrate work here as well.
257
258 Based on odeint from Numerical Recipes (Press, et.al.).
259 If the function you are trying to integrate is not very
260 smooth, or your X values are closely spaced, rk_integrate
261 will probably work better than bs_integrate.
262
263 SEE ALSO: bstoer, rk_integrate, rk_maxits, rk_minstep,
264 rk_maxstep, rk_ngood, rk_nbad, rkdumb, rk4
265 */
266 {
267 if (numberof(rk_x)<2) return array(y1, dimsof(rk_x));
268 local rk_y;
269 if (is_void(dx1)) dx1= rk_x(2)-rk_x(1);
270 rk_nstore= -1;
271 bstoer, derivative, y1, rk_x(1), rk_x(0), epsilon, dx1;
272 return rk_y;
273 }
274
bstoer(derivative,y0,x0,x1,epsilon,dx0)275 func bstoer(derivative, y0,x0, x1,epsilon, dx0)
276 /* DOCUMENT y1= bstoer(derivative, y0,x0, x1,epsilon, dx0)
277 Bulirsch-Stoer integrator, otherwise identical to rkutta routine.
278 All of the options for rkutta (rk_nstore, etc.) work here as well.
279
280 If the function you are trying to integrate is not very
281 smooth, rkutta will probably work better than bstoer.
282
283 SEE ALSO: rkutta, rk_nstore, rk_maxits, rk_minstep, rk_maxstep,
284 rk_ngood, rk_nbad
285 */
286 {
287 extern _rzextr_x, _rzextr_d;
288 rkqc= bsstep;
289 _rzextr_x= array(0.0, numberof(_bs_nseq));
290 _rzextr_d= array(0.+y0, 7);
291 return rkutta(derivative, y0,x0, x1,epsilon, dx0);
292 }
293
294 func bsstep(y,dydx, x,dx, derivative, epsilon,yscale, &dxdid,&dxnxt)
295 {
296 x0= x;
297 y0= y;
298
299 for (;;) {
300 for (i=1 ; i<=numberof(_bs_nseq) ; ++i) {
301 n= _bs_nseq(i);
302 y= rzextr(i, (dx/n)^2, mod_midpt(y0,dydx, x0,dx, derivative, n),
303 yerr, 7);
304 err= max(abs(yerr/yscale))/epsilon;
305 if (err < 1.0) {
306 dxdid= dx;
307 if (i==7) dxnxt= 0.95*dx;
308 else if (i==6) dxnxt= 1.2*dx;
309 else dxnxt= (16./*_bs_nseq(6)*/*dx)/n;
310 return y;
311 }
312 }
313 /* step failed, claimed to be unusual */
314 dx*= 0.0625; /* related to numberof(_bs_nseq) */
315 if (x+dx == x) error, "integration step crash";
316 }
317 }
318
319 _bs_nseq= [2, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96];
320
mod_midpt(y,dydx,x,dx,derivative,nstep)321 func mod_midpt(y,dydx, x,dx, derivative, nstep)
322 {
323 dx/= nstep;
324 ym= y;
325 y+= dydx*dx;
326 x+= dx;
327 dydx= derivative(y,x);
328 dx2= 2.*dx;
329 for (--nstep ; nstep ; --nstep) {
330 swap= ym + dydx*dx2;
331 ym= y;
332 y= swap;
333 x+= dx;
334 dydx= derivative(y,x);
335 }
336 return 0.5*(ym+y+dydx*dx);
337 }
338
339 func rzextr(iest, xest, yest, &yerr, nuse)
340 {
341 _rzextr_x(iest)= xest;
342 if (iest==1) {
343 _rzextr_d(..,1)= yest;
344 yerr= yest;
345 return yest;
346 } else {
347 m1= ((iest<nuse)? iest : nuse) - 1;
348 fx= _rzextr_x(iest-1 : iest-m1 : -1)/xest; /* no more than nuse-1 */
349 yy= yest;
350 v= _rzextr_d(..,1);
351 _rzextr_d(..,1)= c= yy;
352 for (k=1 ; k<=m1 ; ++k) {
353 b1= fx(k)*v;
354 b= b1-c;
355 ok= double(b!=0.0);
356 bad= 1.0-ok;
357 b= ((c-v)/(b+bad))*ok;
358 ddy= c*b + bad*v;
359 c= b1*b + bad*c;
360 if (k!=m1) v= _rzextr_d(..,k+1);
361 _rzextr_d(..,k+1)= ddy;
362 yy+= ddy;
363 }
364 yerr= ddy;
365 return yy;
366 }
367 }
368
369 /* ------------------------------------------------------------------------ */
370