Annotation of OpenXM/src/asir-contrib/packages/src/taka_runge_kutta.rr, Revision 1.23
1.23 ! takayama 1: /* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_runge_kutta.rr,v 1.22 2013/12/15 10:52:36 takayama Exp $ */
1.1 takayama 2: /* From misc/200205/runge-kutta.rr */
3:
4: /* They have not yet been registered in names.rr */
5:
6: #define DEVAL(a) eval(a)
7: Taka_Runge_kutta_adapted0 = 0$
8: Taka_Runge_kutta_epsilon = 0.1$
9: Taka_Runge_kutta_H_Upper_Bound = 0.2$
1.3 takayama 10: Taka_Runge_kutta_Make_Larger = 1$ /* Default 1 */
1.1 takayama 11:
12: Taka_Runge_kutta_graphic0 = 0$ /* load("glib"); */
13: Taka_Runge_kutta_yrange = 10$
14:
15: Taka_Runge_kutta_save_data = 1$
1.15 takayama 16: Taka_Runge_kutta_debug = 0$
1.1 takayama 17:
1.19 takayama 18: extern Tk_rk_deep_eval$
19: Tk_rk_deep_eval=0$
20:
1.1 takayama 21: def taka_runge_kutta_2(F,X,Y,X0,Y0,H,X1) {
22: extern Taka_Runge_kutta_graphic0, Taka_Runge_kutta_yrange, Taka_Runge_kutta_save_data;
23: Alpha =0.5;
24: Beta = 0.5;
25: P = 1; Q = 1;
26:
27: Ans = [];
28: if (Taka_Runge_kutta_graphic0) {
29: glib_open();
30: glib_window(X0,Y0[0]-Taka_Runge_kutta_yrange,X1,Y0[0]+Taka_Runge_kutta_yrange);
31: }
32:
33: if (type(F) == 5) {
34: N = size(F)[0];
35: }else{
36: N = length(F);
37: }
38: if (type(Y0) != 5) {
39: Y0 = newvect(N,Y0);
40: }
41: Yk = Y0;
42: K1 = newvect(N);
43: K2 = newvect(N);
44: Yk1 = newvect(N);
45: Xk = X0;
46:
47: while (Xk < X1) {
48: taka_runge_kutta_replace(K1,F,Y,N,X,Xk,Yk);
49: taka_runge_kutta_replace(K2,F,Y,N,X,Xk+P*H,Yk+Q*H*K1);
50: Yk1 = Yk+H*Alpha*K1+H*Beta*K2;
51: if (Taka_Runge_kutta_save_data) {
52: Ans = cons(cons(Xk,vtol(Yk)),Ans);
53: }
54: print([Xk,Yk[0]]);
55: if (Taka_Runge_kutta_graphic0) glib_line(Xk,Yk[0],Xk+H,Yk1[0]);
56: Xk += H;
57: Yk = Yk1;
58: }
59: return Ans;
60: }
61:
62: def taka_runge_kutta_2_test() {
63: /* Equation of oscilations */
64: F = newvect(2,[y2,-y1]);
65: Y = [y1,y2];
66: taka_runge_kutta_2(F,x,Y,0,[1,0],0.1,15);
67: }
68:
69: def taka_runge_kutta_replace(V,F,Y,N,X,Xk,Rule_vector) {
70: for (I=0; I<N; I++) {
71: V[I] = subst(F[I],X,Xk);
72: for (J=0; J<N; J++) {
73: V[I] = subst(V[I],Y[J],Rule_vector[J]);
74: }
1.9 takayama 75: V[I] = eval(V[I]*exp(0));
1.1 takayama 76: }
77: }
78:
79: def taka_runge_kutta_abs(V) {
80: if (type(V) != 5 && type(V) != 4) { /* not a vector */
81: if (ntype(V) == 4) { /* complex number */
82: return V*conj(V);
83: }else{
84: return(V*V);
85: }
86: }
87: if (type(V) == 5) N = size(V)[0];
88: if (type(V) == 4) N = length(V);
89: S = 0;
90: for (I=0; I<N; I++) {
91: if (ntype(V[I]) == 4) /* complex number */
92: S += V[I]*conj(V[I]);
93: else
94: S += V[I]*V[I];
95: }
96: return S;
97: }
98:
99: def taka_runge_kutta_4(F,X,Y,X0,Y0,H,X1) {
100: /* N : rank of the ODE. */
101: extern Taka_Runge_kutta_adapted0, Taka_Runge_kutta_epsilon,
102: Taka_Runge_kutta_graphic0, Taka_Runge_kutta_yrange,
1.12 takayama 103: Taka_Runge_kutta_save_data, Taka_Runge_kutta_debug;
1.1 takayama 104:
1.10 takayama 105: OneStep = getopt(onestep);
106: if (type(OneStep) <= 0) OneStep = 0; else OneStep = 1;
107: if (OneStep) X1=X0+2*H;
1.14 takayama 108: if ((H<0) && (X1-X0)>0) error("taka_runge_kutta_4, X1-X0 should be <0");
109: if ((H>0) && (X1-X0)<0) error("taka_runge_kutta_4, X1-X0 should be >0");
1.1 takayama 110: Ans = [];
111: if (Taka_Runge_kutta_graphic0) {
112: glib_open();
113: glib_window(X0,Y0[0]-Taka_Runge_kutta_yrange,X1,Y0[0]+Taka_Runge_kutta_yrange);
114: }
1.12 takayama 115: if (X0==X1) return([cons(X0,Y0)]);
1.1 takayama 116:
117: if (type(F) == 5) {
118: N = size(F)[0];
119: }else{
120: N = length(F);
121: }
122: if (type(Y0) != 5) {
123: Y0 = newvect(N,Y0);
124: }
125: Yk = Y0;
126: K1 = newvect(N);
127: K2 = newvect(N);
128: K3 = newvect(N);
129: K4 = newvect(N);
130: Yk1 = newvect(N);
131: Xk = X0;
132:
1.14 takayama 133: while (H<0? Xk > X1: Xk < X1) {
1.1 takayama 134: taka_runge_kutta_replace(K1,F,Y,N,X,Xk,Yk);
1.3 takayama 135: taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H*(1/2),Yk+K1*(1/2)*H);
136: taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H*(1/2),Yk+K2*(1/2)*H);
1.1 takayama 137: taka_runge_kutta_replace(K4,F,Y,N,X,Xk+H,Yk+K3*H);
138: Yk1 = Yk+H*(K1/6+K2/3+K3/3+K4/6);
1.12 takayama 139: if (Taka_Runge_kutta_debug) print([Xk,Yk[0]]);
1.1 takayama 140: if (Taka_Runge_kutta_save_data) {
141: Ans = cons(cons(Xk,vtol(Yk)),Ans);
142: }
1.10 takayama 143: if (OneStep) {
1.12 takayama 144: return([cons(Xk+H,vtol(Yk1)), cons(Xk,vtol(Yk))]);
1.10 takayama 145: }
1.1 takayama 146: if (Taka_Runge_kutta_graphic0) glib_line(Xk,Yk[0],Xk+H,Yk1[0]);
147: if (Taka_Runge_kutta_adapted0 &&
148: (taka_runge_kutta_abs(Yk1-Yk) > Taka_Runge_kutta_epsilon)) {
1.3 takayama 149: H = H*(1/2);
1.1 takayama 150: }else{
151: if (Taka_Runge_kutta_adapted0) H = H*2;
152: Xk += H;
153: Yk = Yk1;
154: }
155: }
156: return Ans;
157: }
158:
159: def taka_runge_kutta_4_adaptive(F,X,Y,X0,Y0,H,X1) {
160: /* N : rank of the ODE. */
161: extern Taka_Runge_kutta_epsilon,
162: Taka_Runge_kutta_graphic0, Taka_Runge_kutta_yrange,
163: Taka_Runge_kutta_save_data,
164: Taka_Runge_kutta_H_Upper_Bound,
165: Taka_Runge_kutta_Make_Larger;
166:
1.5 takayama 167: Opt = getopt();
1.6 takayama 168: if (taka_runge_kutta_complex_gt(H,0)) Forward = 1; else Forward = 0;
1.5 takayama 169: while(Opt != []) {
170: if (car(Opt)[0] == "forward") {
171: Forward = car(Opt)[1];
172: }
173: Opt = cdr(Opt);
174: }
175:
1.1 takayama 176: Ans = [cons(X0,Y0)];
177: if (Taka_Runge_kutta_graphic0) {
178: glib_open();
179: glib_window(X0,Y0[0]-Taka_Runge_kutta_yrange,X1,Y0[0]+Taka_Runge_kutta_yrange);
180: }
181:
182: if (type(F) == 5) {
183: N = size(F)[0];
184: }else{
185: N = length(F);
186: }
187: if (type(Y0) != 5) {
188: Y0 = newvect(N,Y0);
189: }
190: Yk = Y0;
191: K1 = newvect(N);
192: K2 = newvect(N);
193: K3 = newvect(N);
194: K4 = newvect(N);
195: Yk1 = newvect(N);
1.3 takayama 196: Yk2 = newvect(N);
197: Yk3 = newvect(N);
1.1 takayama 198: Xk = X0;
199:
200: while (true) {
1.5 takayama 201: if (Forward) {
1.6 takayama 202: /* if (Xk > X1) break; */
203: if (taka_runge_kutta_complex_gt(Xk,X1)) break;
1.5 takayama 204: } else{
1.6 takayama 205: /* if (Xk < X1) break; */
206: if (taka_runge_kutta_complex_gt(X1,Xk)) break;
1.1 takayama 207: }
208: /* Regular step */
209: taka_runge_kutta_replace(K1,F,Y,N,X,Xk,Yk);
1.3 takayama 210: taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H*(1/2),Yk+K1*(1/2)*H);
211: taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H*(1/2),Yk+K2*(1/2)*H);
1.1 takayama 212: taka_runge_kutta_replace(K4,F,Y,N,X,Xk+H,Yk+K3*H);
213: Yk1 = Yk+H*(K1/6+K2/3+K3/3+K4/6);
214: /* half step */
215: H2 = H/2;
216: taka_runge_kutta_replace(K1,F,Y,N,X,Xk,Yk);
1.3 takayama 217: taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H2*(1/2),Yk+K1*(1/2)*H2);
218: taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H2*(1/2),Yk+K2*(1/2)*H2);
219: taka_runge_kutta_replace(K4,F,Y,N,X,Xk+H2,Yk+K3*H2);
1.1 takayama 220: Yk2 = Yk+H2*(K1/6+K2/3+K3/3+K4/6);
221:
1.3 takayama 222: taka_runge_kutta_replace(K1,F,Y,N,X,Xk+H2,Yk2);
223: taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H2+H2*(1/2),Yk2+K1*(1/2)*H2);
224: taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H2+H2*(1/2),Yk2+K2*(1/2)*H2);
225: taka_runge_kutta_replace(K4,F,Y,N,X,Xk+H2+H2,Yk2+K3*H2);
226: Yk3 = Yk2+H2*(K1/6+K2/3+K3/3+K4/6);
227:
228: /* This is a strategy which you may change. */
1.4 takayama 229: /* WantedPrec = Taka_Runge_kutta_epsilon*taka_runge_kutta_abs(Yk);*/
230: WantedPrec = Taka_Runge_kutta_epsilon;
1.3 takayama 231:
232: Delta1 = DEVAL(taka_runge_kutta_abs(Yk3-Yk1));
233: if (Delta1 != 0) {
234: Habs = DEVAL((WantedPrec/Delta1)^(1/5));
235: Habs = (4/5)*Habs; /* 0.8 = (4/5) is the safety factor */
236: }else{
237: Habs = 2; /* Any large number */
238: }
1.4 takayama 239: /* print("Habs="+rtostr(Habs)); */
1.3 takayama 240: if (Habs < 1) { /* Compute again. */
241: H = H*Habs;
1.1 takayama 242: print("Changing to Smaller step size: "+rtostr(H));
243: print([Xk,Yk[0]]);
244: }else{ /* Go ahead */
245: Xk += H;
246: Yk = Yk1;
1.3 takayama 247: if ((H<Taka_Runge_kutta_H_Upper_Bound) && Taka_Runge_kutta_Make_Larger) {
248: H = (Habs*H > Taka_Runge_kutta_H_Upper_Bound?
1.4 takayama 249: (H/number_abs(H))*Taka_Runge_kutta_H_Upper_Bound :
1.3 takayama 250: Habs*H); /* Habs*H2*2 */
1.1 takayama 251: print("Changing to a larger step size: "+rtostr(H));
252: }
253: print([Xk,Yk[0]]);
254: if (Taka_Runge_kutta_save_data) {
255: Ans = cons(cons(Xk,vtol(Yk)),Ans);
256: }
257: if (Taka_Runge_kutta_graphic0) glib_line(Xk,Yk[0],Xk+H,Yk1[0]);
258: }
259: }
260: return Ans;
261: }
262:
263: /* load("glib"); load("taka_plot.rr"); to execute the functions below. */
264: def taka_runge_kutta_4_a_test() {
265: /* exponential function */
266: F = newvect(1,[y1]);
267: Y = [y1];
1.3 takayama 268: T = taka_runge_kutta_4_adaptive(F,x,Y,0,[1],0.1,5);
269: taka_plot_auto(T);
1.13 takayama 270: print("Eval by eval(exp(?)) : ",0); print([T[0][0],eval(exp(T[0][0]))]);
1.3 takayama 271: }
272:
273: def taka_runge_kutta_4_a2_test() {
274: /* Equation of oscilations */
275: F = newvect(2,[y2,-y1]);
276: Y = [y1,y2];
277: T = taka_runge_kutta_4_adaptive(F,x,Y,0,[1,0],0.1,15);
1.1 takayama 278: taka_plot_auto(T);
1.13 takayama 279: print("Eval by eval((?)) : ",0); print([T[0][0],eval(cos(T[0][0]))]);
1.1 takayama 280: }
281:
282: def taka_runge_kutta_4_test() {
283: /* Equation of oscilations */
284: F = newvect(2,[y2,-y1]);
285: Y = [y1,y2];
286: T=taka_runge_kutta_4(F,x,Y,0,[1,0],0.1,15);
287: print(T);
288: taka_plot_auto(T);
289: }
290:
1.14 takayama 291: def taka_runge_kutta_4_test2() {
292: /* Equation of oscilations */
293: F = newvect(2,[y2,-y1]);
294: Y = [y1,y2];
295: T=taka_runge_kutta_4(F,x,Y,15,[1,0],-0.1,0);
296: print(T);
297: taka_plot_auto(T);
298: }
299:
1.21 takayama 300: def taka_runge_kutta_replace_linear(V,F,Y,N,X,Xk,Rule_vector,Inhom) {
1.19 takayama 301: extern Tk_rk_deep_eval;
302: if (Tk_rk_deep_eval) {
303: V1=base_replace(F,[[X,Xk]]);
304: }else{
305: V1=base_replace_n(F,[[X,Xk]]);
306: }
1.15 takayama 307: V1=V1*Rule_vector;
308: for (I=0; I<N; I++) {
1.19 takayama 309: if (Tk_rk_deep_eval)
310: V[I] = eval(V1[I]*exp(0));
311: else
312: V[I] = V1[I];
1.15 takayama 313: }
1.21 takayama 314: if (Inhom != 0) {
315: /* printf("X=%a, Xk=%a\n",X,Xk); */
316: if (Tk_rk_deep_eval) {
317: V1=map(eval,base_replace(Inhom,[[X,Xk]]));
318: }else{
319: V1=base_replace_n(Inhom,[[X,Xk]]);
320: }
321: for (I=0; I<N; I++) V[I] += V1[I];
322: }
1.15 takayama 323: }
324:
325: /* Y is a dummy */
326: def taka_runge_kutta_4_linear(F,X,Y,X0,Y0,H,X1) {
327: /* N : rank of the ODE. */
328: extern Taka_Runge_kutta_adapted0, Taka_Runge_kutta_epsilon,
329: Taka_Runge_kutta_graphic0, Taka_Runge_kutta_yrange,
330: Taka_Runge_kutta_save_data, Taka_Runge_kutta_debug;
331:
332: OneStep = getopt(onestep);
333: if (type(OneStep) <= 0) OneStep = 0; else OneStep = 1;
334: if (OneStep) X1=X0+2*H;
1.21 takayama 335: Inhom = getopt(inhom);
336: if (type(Inhom) < 0) Inhom = 0;
337: /* else printf("inhom=%a, onestep=%a\n",Inhom,OneStep); */ /* for debug */
338: if (type(Inhom) == 4) {
339: Inhom = newvect(length(Inhom),Inhom);
340: }
1.15 takayama 341: if ((H<0) && (X1-X0)>0) error("taka_runge_kutta_4_linear, X1-X0 should be <0");
342: if ((H>0) && (X1-X0)<0) error("taka_runge_kutta_4_linear, X1-X0 should be >0");
343: Ans = [];
344: if (Taka_Runge_kutta_graphic0) {
345: glib_open();
346: glib_window(X0,Y0[0]-Taka_Runge_kutta_yrange,X1,Y0[0]+Taka_Runge_kutta_yrange);
347: }
348: if (X0==X1) return([cons(X0,Y0)]);
349:
350: if (type(F) == 4) {
351: F=newmat(length(F),length(F[0]),F);
352: }
353: N = size(F)[0];
354:
355: if (type(Y0) != 5) {
356: Y0 = newvect(N,Y0);
357: }
358: Yk = Y0;
359: K1 = newvect(N);
360: K2 = newvect(N);
361: K3 = newvect(N);
362: K4 = newvect(N);
363: Yk1 = newvect(N);
364: Xk = X0;
365:
366: while (H<0? Xk > X1: Xk < X1) {
1.21 takayama 367: taka_runge_kutta_replace_linear(K1,F,Y,N,X,Xk,Yk,Inhom);
368: taka_runge_kutta_replace_linear(K2,F,Y,N,X,Xk+H*(1/2),Yk+K1*(1/2)*H,Inhom);
369: taka_runge_kutta_replace_linear(K3,F,Y,N,X,Xk+H*(1/2),Yk+K2*(1/2)*H,Inhom);
370: taka_runge_kutta_replace_linear(K4,F,Y,N,X,Xk+H,Yk+K3*H,Inhom);
1.15 takayama 371: Yk1 = Yk+H*(K1/6+K2/3+K3/3+K4/6);
372: if (Taka_Runge_kutta_debug) print([Xk,Yk[0]]);
373: if (Taka_Runge_kutta_save_data) {
374: Ans = cons(cons(Xk,vtol(Yk)),Ans);
375: }
376: if (OneStep) {
377: return([cons(Xk+H,vtol(Yk1)), cons(Xk,vtol(Yk))]);
378: }
379: if (Taka_Runge_kutta_graphic0) glib_line(Xk,Yk[0],Xk+H,Yk1[0]);
380: if (Taka_Runge_kutta_adapted0 &&
381: (taka_runge_kutta_abs(Yk1-Yk) > Taka_Runge_kutta_epsilon)) {
382: H = H*(1/2);
383: }else{
384: if (Taka_Runge_kutta_adapted0) H = H*2;
385: Xk += H;
386: Yk = Yk1;
387: }
388: }
389: return Ans;
390: }
391:
392: def taka_runge_kutta_4_linear_test() {
393: /* Airy equation y''-x y = 0
394: [evalf(AiryAi(0)),evalf(subs(x=0,diff(AiryAi(x),x)))];
395: Y0=[0.3550280540, -0.2588194038]
396: evalf(AiryAi(-5)); --> 0.35076
397: */
398: F = [[0,1],[x,0]];
399: Y = [y1,y2];
400: Y0=[0.3550280540, -0.2588194038];
401: T=taka_runge_kutta_4_linear(F,x,Y,0,Y0,-0.1,-5);
402: print(T);
403: taka_plot_auto(T);
404: T2=taka_runge_kutta_4([y2,x*y1],x,Y,0,Y0,-0.1,-5);
405: print("AiryAi(-5) --> 0.35076");
406: return([T[0],T2[0]]);
407: }
408:
409: /*
410: def base_replace_n(F,R) { return base_replace(F,R); }
411: */
1.6 takayama 412:
413: /* cf. asir2000/engine/cplx.c int cmpcplx(a,b),
414: which does not compare the real part and imaginary part.
415: Instead, it compares NID (number id)
416: */
417: def taka_runge_kutta_complex_gt(A,B) {
418: Ar = number_real_part(A); Ai = number_imaginary_part(A);
419: Br = number_real_part(B); Bi = number_imaginary_part(B);
420: if (Ar > Br) return 1;
421: else if (Ar < Br) return 0;
422: if (Ai > Bi) return 1;
423: else if (Ai < Bi) return 0;
424: return 0;
425: }
1.1 takayama 426:
427: Loaded_taka_runge_kutta=1$
1.7 takayama 428:
429: /* cf. misc-2003/09/neval/ellip.* */
430: /* runge_kutta_4 is still buggy for complex numbers */
431:
1.16 takayama 432: module tk_rk;
433: localf taka_minus;
434: localf taka_runge_kutta_reverse ;
435: localf taka_runge_kutta_4a;
436: localf taka_runge_kutta_4a_linear;
437: localf test4 ;
438: localf test4b ;
1.18 takayama 439: localf runge_kutta_4;
440: localf runge_kutta_4_linear;
1.16 takayama 441: def taka_minus(Ob) {
1.17 takayama 442: if (type(Ob) != 4) return(-Ob);
443: else return map(taka_minus,Ob);
1.16 takayama 444: }
445: def taka_runge_kutta_reverse(A) {
446: B=[];
447: for (; length(A) != 0; A=cdr(A)) {
448: T=car(A);
449: B=cons(cons(-T[0],cdr(T)),B);
450: }
451: return reverse(B);
452: }
453: def taka_runge_kutta_4a(FF,X0,Y,S0,Ys,T0,H) {
454: if (T0 < S0) {
455: /* opposite direction */
456: return taka_runge_kutta_reverse(
457: taka_runge_kutta_4a(map(taka_minus,base_replace(FF,[[X0,-X0]])),X0,Y,-S0,Ys,-T0,H));
458: }
459: if (H >= T0-S0) {
460: A=taka_runge_kutta_4(FF,X0,Y,S0,Ys,T0-S0,0 | onestep=1);
461: }else{
462: A=taka_runge_kutta_4(FF,X0,Y,S0,Ys,H,T0);
463: T=A[0];
464: if (T0-T[0] > 0) {
465: B=taka_runge_kutta_4(FF,X0,Y,T[0],cdr(T),T0-T[0],0 | onestep=1);
466: T=B[0];
467: A=cons(T,A);
468: }
469: }
470: return(A);
471: }
1.18 takayama 472: def runge_kutta_4(FF,X0,Y,S0,Ys,T0,H) {
473: return taka_runge_kutta_4a(FF,X0,Y,S0,Ys,T0,H);
474: }
475: def runge_kutta_4_linear(FF,X0,Y,S0,Ys,T0,H) {
1.21 takayama 476: Opt = getopt();
477: return taka_runge_kutta_4a_linear(FF,X0,Y,S0,Ys,T0,H | option_list=Opt);
1.18 takayama 478: }
1.16 takayama 479:
480: def taka_runge_kutta_4a_linear(FF,X0,Y,S0,Ys,T0,H) {
1.21 takayama 481: Opt = getopt();
1.17 takayama 482: if (T0 < S0) {
483: /* opposite direction */
1.20 takayama 484: if (H<0) H=-H;
1.22 takayama 485: Opt = [];
486: if (getopt(onestep) >= 0) Opt=[["onestep",getopt(onestep)]];
487: if (getopt(inhom) >= 0) Inhom = getopt(inhom);
488: if (Inhom != 0) {
489: Inhom = map(taka_minus,base_replace(Inhom,[[X0,-X0]]));
490: Opt = append(Opt,[["inhom",Inhom]]);
491: }
1.17 takayama 492: return taka_runge_kutta_reverse(
1.21 takayama 493: taka_runge_kutta_4a_linear(map(taka_minus,base_replace(FF,[[X0,-X0]])),X0,Y,-S0,Ys,-T0,H | option_list=Opt));
1.17 takayama 494: }
1.16 takayama 495: if (H >= T0-S0) {
1.21 takayama 496: A=taka_runge_kutta_4_linear(FF,X0,Y,S0,Ys,T0-S0,0 | option_list=append(Opt,[["onestep",1]]));
1.16 takayama 497: }else{
1.21 takayama 498: A=taka_runge_kutta_4_linear(FF,X0,Y,S0,Ys,H,T0 | option_list=Opt);
1.16 takayama 499: T=A[0];
500: if (T0-T[0] > 0) {
1.23 ! takayama 501: B=taka_runge_kutta_4_linear(FF,X0,Y,T[0],cdr(T),T0-T[0],T0+(T0-T[0])/10 | option_list=append(Opt,[["onestep",1]]));
1.16 takayama 502: T=B[0];
503: A=cons(T,A);
504: }
505: }
506: return(A);
507: }
508:
509: /* equation of oscilation */
510: def test4() {
1.18 takayama 511: A=runge_kutta_4([y1,-y0],x,[y0,y1],0,[1,0],3.14*2,0.1);
1.16 takayama 512: taka_plot_auto(A);
513: return(A);
514: }
515:
516: def test4b() {
517: A=taka_runge_kutta_4a([y1,-y0],x,[y0,y1],3.14,[-1,0],0,0.1);
518: taka_plot_auto(A);
519: return(A);
520: }
521: endmodule;
522:
1.7 takayama 523: import("taka_plot.rr")$
524: pari(allocatemem,10^7)$
525: module rktest;
526: localf re$
527: localf re2$
528: localf im$
529: localf im2$
530: localf tryA$
531: localf tryA2$
532: localf geq$
533: localf test1$
1.21 takayama 534: localf rcheck1$
1.7 takayama 535:
536: def re(L) {
537: return map(re2,L);
538: }
539: def re2(P) {
540: return map(number_real_part,P);
541: }
542: def im(L) {
543: return map(im2,L);
544: }
545: def im2(P) {
546: return map(number_imaginary_part,P);
547: }
548:
549: def geq() {
1.8 takayama 550: L=x*(1-x)*dx^2+(c-(a+b+1)*x)*dx-a*b;
1.7 takayama 551: L=base_replace(L,[[a,1/2],[b,1/2],[c,1]]);
552:
1.8 takayama 553: L2 = -((c-(a+b+1)*x)*y2-a*b*y1)/(x*(1-x));
1.7 takayama 554: L2=base_replace(L2,[[a,1/2],[b,1/2],[c,1]]);
555: return [ y2, L2];
556: }
557: def tryA() {
558: LL = geq();
559: A = taka_runge_kutta_4_adaptive(
560: LL,
561: x,[y1,y2],
562: 0.5+0.5*@i,[-6.78383-1.28991*@i, -1.51159-1.7935*@i],
563: (3-@i)*0.0005, 2.0);
564: taka_plot_auto(re(A));
565: return A;
566: }
567: def tryA2() {
568: LL = geq();
569: A = taka_runge_kutta_4(
570: LL,
571: x,[y1,y2],
572: 0.5+0.5*@i,[-6.78383-1.28991*@i, -1.51159-1.7935*@i],
573: (3-@i)*0.0005, 1.0);
574: taka_plot_auto(re(A));
575: return A;
576: }
577:
578: def test1() {
579: A=tryA2();
580: B=A[100];
581: print(B);
1.10 takayama 582: /* cf. misc-2008/A2/misc/ellip2.m
1.7 takayama 583: p2 = 0.5+0.5*I --> [-6.78383-1.28991*@i, -1.51159-1.7935*@i].
1.8 takayama 584: p2 = 0.8495+0.3835*I;
1.7 takayama 585: Print["-------------------"];
586: Print[p2];
587: Print[N[-2*Gamma[1/2]^2*Hypergeometric2F1[1/2,1/2,1,z] /. {z->p2}]]
588: Print[N[D[-2*Gamma[1/2]^2*Hypergeometric2F1[1/2,1/2,1,z],z] /. {z->p2}]]
589: */
1.8 takayama 590: print("math: [0.8495+0.3835*I, -7.64079 - 2.0799*I, -1.16364 - 4.14709*I] ");
1.10 takayama 591: print("It was Buggy tryA2() and tryA() --> fixed. see log of 1.8");
1.7 takayama 592: }
593:
1.21 takayama 594: def rcheck1() {
595: /* f'=f/z+z*cos(z) where f=z*sin(z) */
596: Tk_rk_deep_eval=1;
597: Step = 0.01;
598: X0 = 3.14/2;
599: X = 3.14;
600: Iv=[X0*1];
601: Eq=[[1/z]]; Inhom = [z*cos(z)];
602: A=tk_rk.runge_kutta_4_linear(Eq,z,[],X0,Iv,X,Step | inhom=Inhom, onestep=1);
603: print(A);
604: A=tk_rk.runge_kutta_4_linear(Eq,z,[],X0,Iv,X,Step | inhom=Inhom);
605: return A;
606: }
607:
1.7 takayama 608: endmodule;
609:
1.1 takayama 610: end$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to taka_runge_kutta.rr CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / packages / src