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