=================================================================== RCS file: /home/cvs/OpenXM/src/asir-contrib/packages/src/taka_runge_kutta.rr,v retrieving revision 1.1 retrieving revision 1.18 diff -u -p -r1.1 -r1.18 --- OpenXM/src/asir-contrib/packages/src/taka_runge_kutta.rr 2002/11/16 08:09:31 1.1 +++ OpenXM/src/asir-contrib/packages/src/taka_runge_kutta.rr 2010/09/29 00:33:24 1.18 @@ -1,4 +1,4 @@ -/* $Id: taka_runge_kutta.rr,v 1.1 2002/11/16 08:09:31 takayama Exp $ */ +/* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_runge_kutta.rr,v 1.17 2010/04/18 01:08:37 takayama Exp $ */ /* From misc/200205/runge-kutta.rr */ /* They have not yet been registered in names.rr */ @@ -7,14 +7,14 @@ Taka_Runge_kutta_adapted0 = 0$ Taka_Runge_kutta_epsilon = 0.1$ Taka_Runge_kutta_H_Upper_Bound = 0.2$ -Taka_Runge_kutta_Make_Larger = 0$ /* Default 1, H = 2*H */ +Taka_Runge_kutta_Make_Larger = 1$ /* Default 1 */ Taka_Runge_kutta_graphic0 = 0$ /* load("glib"); */ Taka_Runge_kutta_yrange = 10$ Taka_Runge_kutta_save_data = 1$ +Taka_Runge_kutta_debug = 0$ - def taka_runge_kutta_2(F,X,Y,X0,Y0,H,X1) { extern Taka_Runge_kutta_graphic0, Taka_Runge_kutta_yrange, Taka_Runge_kutta_save_data; Alpha =0.5; @@ -69,6 +69,7 @@ def taka_runge_kutta_replace(V,F,Y,N,X,Xk,Rule_vector) for (J=0; J0) error("taka_runge_kutta_4, X1-X0 should be <0"); + if ((H>0) && (X1-X0)<0) error("taka_runge_kutta_4, X1-X0 should be >0"); Ans = []; if (Taka_Runge_kutta_graphic0) { glib_open(); glib_window(X0,Y0[0]-Taka_Runge_kutta_yrange,X1,Y0[0]+Taka_Runge_kutta_yrange); } + if (X0==X1) return([cons(X0,Y0)]); if (type(F) == 5) { N = size(F)[0]; @@ -120,20 +127,23 @@ def taka_runge_kutta_4(F,X,Y,X0,Y0,H,X1) { Yk1 = newvect(N); Xk = X0; - while (Xk < X1) { + while (H<0? Xk > X1: Xk < X1) { taka_runge_kutta_replace(K1,F,Y,N,X,Xk,Yk); - taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H*0.5,Yk+K1*0.5*H); - taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H*0.5,Yk+K2*0.5*H); + taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H*(1/2),Yk+K1*(1/2)*H); + taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H*(1/2),Yk+K2*(1/2)*H); taka_runge_kutta_replace(K4,F,Y,N,X,Xk+H,Yk+K3*H); Yk1 = Yk+H*(K1/6+K2/3+K3/3+K4/6); - print([Xk,Yk[0]]); + if (Taka_Runge_kutta_debug) print([Xk,Yk[0]]); if (Taka_Runge_kutta_save_data) { Ans = cons(cons(Xk,vtol(Yk)),Ans); } + if (OneStep) { + return([cons(Xk+H,vtol(Yk1)), cons(Xk,vtol(Yk))]); + } if (Taka_Runge_kutta_graphic0) glib_line(Xk,Yk[0],Xk+H,Yk1[0]); if (Taka_Runge_kutta_adapted0 && (taka_runge_kutta_abs(Yk1-Yk) > Taka_Runge_kutta_epsilon)) { - H = H*0.5; + H = H*(1/2); }else{ if (Taka_Runge_kutta_adapted0) H = H*2; Xk += H; @@ -151,6 +161,15 @@ def taka_runge_kutta_4_adaptive(F,X,Y,X0,Y0,H,X1) { Taka_Runge_kutta_H_Upper_Bound, Taka_Runge_kutta_Make_Larger; + Opt = getopt(); + if (taka_runge_kutta_complex_gt(H,0)) Forward = 1; else Forward = 0; + while(Opt != []) { + if (car(Opt)[0] == "forward") { + Forward = car(Opt)[1]; + } + Opt = cdr(Opt); + } + Ans = [cons(X0,Y0)]; if (Taka_Runge_kutta_graphic0) { glib_open(); @@ -171,42 +190,61 @@ def taka_runge_kutta_4_adaptive(F,X,Y,X0,Y0,H,X1) { K3 = newvect(N); K4 = newvect(N); Yk1 = newvect(N); + Yk2 = newvect(N); + Yk3 = newvect(N); Xk = X0; while (true) { - if (H > 0) { - if (Xk > X1) break; - }else{ - if (Xk < X1) break; + if (Forward) { + /* if (Xk > X1) break; */ + if (taka_runge_kutta_complex_gt(Xk,X1)) break; + } else{ + /* if (Xk < X1) break; */ + if (taka_runge_kutta_complex_gt(X1,Xk)) break; } /* Regular step */ taka_runge_kutta_replace(K1,F,Y,N,X,Xk,Yk); - taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H*0.5,Yk+K1*0.5*H); - taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H*0.5,Yk+K2*0.5*H); + taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H*(1/2),Yk+K1*(1/2)*H); + taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H*(1/2),Yk+K2*(1/2)*H); taka_runge_kutta_replace(K4,F,Y,N,X,Xk+H,Yk+K3*H); Yk1 = Yk+H*(K1/6+K2/3+K3/3+K4/6); /* half step */ H2 = H/2; taka_runge_kutta_replace(K1,F,Y,N,X,Xk,Yk); - taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H2*0.5,Yk+K1*0.5*H2); - taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H2*0.5,Yk+K2*0.5*H2); - taka_runge_kutta_replace(K4,F,Y,N,X,Xk+H,Yk+K3*H2); + taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H2*(1/2),Yk+K1*(1/2)*H2); + taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H2*(1/2),Yk+K2*(1/2)*H2); + taka_runge_kutta_replace(K4,F,Y,N,X,Xk+H2,Yk+K3*H2); Yk2 = Yk+H2*(K1/6+K2/3+K3/3+K4/6); - Delta1 = DEVAL(taka_runge_kutta_abs(Yk2-Yk1))^(1/2); - Habs = DEVAL(taka_runge_kutta_abs(H)^(1/2)); - HHabs = DEVAL(Habs*exp(0.2*log(Taka_Runge_kutta_epsilon/Delta1))); - /* print(HH); */ - if (HHabs < Habs) { /* Compute again. */ - H = H/2; + taka_runge_kutta_replace(K1,F,Y,N,X,Xk+H2,Yk2); + taka_runge_kutta_replace(K2,F,Y,N,X,Xk+H2+H2*(1/2),Yk2+K1*(1/2)*H2); + taka_runge_kutta_replace(K3,F,Y,N,X,Xk+H2+H2*(1/2),Yk2+K2*(1/2)*H2); + taka_runge_kutta_replace(K4,F,Y,N,X,Xk+H2+H2,Yk2+K3*H2); + Yk3 = Yk2+H2*(K1/6+K2/3+K3/3+K4/6); + + /* This is a strategy which you may change. */ + /* WantedPrec = Taka_Runge_kutta_epsilon*taka_runge_kutta_abs(Yk);*/ + WantedPrec = Taka_Runge_kutta_epsilon; + + Delta1 = DEVAL(taka_runge_kutta_abs(Yk3-Yk1)); + if (Delta1 != 0) { + Habs = DEVAL((WantedPrec/Delta1)^(1/5)); + Habs = (4/5)*Habs; /* 0.8 = (4/5) is the safety factor */ + }else{ + Habs = 2; /* Any large number */ + } + /* print("Habs="+rtostr(Habs)); */ + if (Habs < 1) { /* Compute again. */ + H = H*Habs; print("Changing to Smaller step size: "+rtostr(H)); print([Xk,Yk[0]]); }else{ /* Go ahead */ Xk += H; Yk = Yk1; - if ((HHabs > 2*Habs) && (H Taka_Runge_kutta_H_Upper_Bound? + (H/number_abs(H))*Taka_Runge_kutta_H_Upper_Bound : + Habs*H); /* Habs*H2*2 */ print("Changing to a larger step size: "+rtostr(H)); } print([Xk,Yk[0]]); @@ -221,18 +259,23 @@ def taka_runge_kutta_4_adaptive(F,X,Y,X0,Y0,H,X1) { /* load("glib"); load("taka_plot.rr"); to execute the functions below. */ def taka_runge_kutta_4_a_test() { - /* Equation of oscilations */ - /* F = newvect(2,[y2,-y1]); - Y = [y1,y2]; - T = taka_runge_kutta_4_adaptive(F,x,Y,0,[1,0],0.1,15); */ /* exponential function */ F = newvect(1,[y1]); Y = [y1]; - T = taka_runge_kutta_4_adaptive(F,x,Y,0,[1],0.1,2); + T = taka_runge_kutta_4_adaptive(F,x,Y,0,[1],0.1,5); taka_plot_auto(T); - print([T[0][0],eval(exp(T[0][0]))]); + print("Eval by eval(exp(?)) : ",0); print([T[0][0],eval(exp(T[0][0]))]); } +def taka_runge_kutta_4_a2_test() { + /* Equation of oscilations */ + F = newvect(2,[y2,-y1]); + Y = [y1,y2]; + T = taka_runge_kutta_4_adaptive(F,x,Y,0,[1,0],0.1,15); + taka_plot_auto(T); + print("Eval by eval((?)) : ",0); print([T[0][0],eval(cos(T[0][0]))]); +} + def taka_runge_kutta_4_test() { /* Equation of oscilations */ F = newvect(2,[y2,-y1]); @@ -242,6 +285,275 @@ def taka_runge_kutta_4_test() { taka_plot_auto(T); } +def taka_runge_kutta_4_test2() { + /* Equation of oscilations */ + F = newvect(2,[y2,-y1]); + Y = [y1,y2]; + T=taka_runge_kutta_4(F,x,Y,15,[1,0],-0.1,0); + print(T); + taka_plot_auto(T); +} +def taka_runge_kutta_replace_linear(V,F,Y,N,X,Xk,Rule_vector) { + V1=base_replace_n(F,[[X,Xk]]); + V1=V1*Rule_vector; + for (I=0; I0) error("taka_runge_kutta_4_linear, X1-X0 should be <0"); + if ((H>0) && (X1-X0)<0) error("taka_runge_kutta_4_linear, X1-X0 should be >0"); + Ans = []; + if (Taka_Runge_kutta_graphic0) { + glib_open(); + glib_window(X0,Y0[0]-Taka_Runge_kutta_yrange,X1,Y0[0]+Taka_Runge_kutta_yrange); + } + if (X0==X1) return([cons(X0,Y0)]); + + if (type(F) == 4) { + F=newmat(length(F),length(F[0]),F); + } + N = size(F)[0]; + + if (type(Y0) != 5) { + Y0 = newvect(N,Y0); + } + Yk = Y0; + K1 = newvect(N); + K2 = newvect(N); + K3 = newvect(N); + K4 = newvect(N); + Yk1 = newvect(N); + Xk = X0; + + while (H<0? Xk > X1: Xk < X1) { + taka_runge_kutta_replace_linear(K1,F,Y,N,X,Xk,Yk); + taka_runge_kutta_replace_linear(K2,F,Y,N,X,Xk+H*(1/2),Yk+K1*(1/2)*H); + taka_runge_kutta_replace_linear(K3,F,Y,N,X,Xk+H*(1/2),Yk+K2*(1/2)*H); + taka_runge_kutta_replace_linear(K4,F,Y,N,X,Xk+H,Yk+K3*H); + Yk1 = Yk+H*(K1/6+K2/3+K3/3+K4/6); + if (Taka_Runge_kutta_debug) print([Xk,Yk[0]]); + if (Taka_Runge_kutta_save_data) { + Ans = cons(cons(Xk,vtol(Yk)),Ans); + } + if (OneStep) { + return([cons(Xk+H,vtol(Yk1)), cons(Xk,vtol(Yk))]); + } + if (Taka_Runge_kutta_graphic0) glib_line(Xk,Yk[0],Xk+H,Yk1[0]); + if (Taka_Runge_kutta_adapted0 && + (taka_runge_kutta_abs(Yk1-Yk) > Taka_Runge_kutta_epsilon)) { + H = H*(1/2); + }else{ + if (Taka_Runge_kutta_adapted0) H = H*2; + Xk += H; + Yk = Yk1; + } + } + return Ans; +} + +def taka_runge_kutta_4_linear_test() { + /* Airy equation y''-x y = 0 + [evalf(AiryAi(0)),evalf(subs(x=0,diff(AiryAi(x),x)))]; + Y0=[0.3550280540, -0.2588194038] + evalf(AiryAi(-5)); --> 0.35076 + */ + F = [[0,1],[x,0]]; + Y = [y1,y2]; + Y0=[0.3550280540, -0.2588194038]; + T=taka_runge_kutta_4_linear(F,x,Y,0,Y0,-0.1,-5); + print(T); + taka_plot_auto(T); + T2=taka_runge_kutta_4([y2,x*y1],x,Y,0,Y0,-0.1,-5); + print("AiryAi(-5) --> 0.35076"); + return([T[0],T2[0]]); +} + +/* +def base_replace_n(F,R) { return base_replace(F,R); } +*/ + +/* cf. asir2000/engine/cplx.c int cmpcplx(a,b), + which does not compare the real part and imaginary part. + Instead, it compares NID (number id) +*/ +def taka_runge_kutta_complex_gt(A,B) { + Ar = number_real_part(A); Ai = number_imaginary_part(A); + Br = number_real_part(B); Bi = number_imaginary_part(B); + if (Ar > Br) return 1; + else if (Ar < Br) return 0; + if (Ai > Bi) return 1; + else if (Ai < Bi) return 0; + return 0; +} + Loaded_taka_runge_kutta=1$ + +/* cf. misc-2003/09/neval/ellip.* */ +/* runge_kutta_4 is still buggy for complex numbers */ + +module tk_rk; +localf taka_minus; +localf taka_runge_kutta_reverse ; +localf taka_runge_kutta_4a; +localf taka_runge_kutta_4a_linear; +localf test4 ; +localf test4b ; +localf runge_kutta_4; +localf runge_kutta_4_linear; +def taka_minus(Ob) { + if (type(Ob) != 4) return(-Ob); + else return map(taka_minus,Ob); +} +def taka_runge_kutta_reverse(A) { + B=[]; + for (; length(A) != 0; A=cdr(A)) { + T=car(A); + B=cons(cons(-T[0],cdr(T)),B); + } + return reverse(B); +} +def taka_runge_kutta_4a(FF,X0,Y,S0,Ys,T0,H) { + if (T0 < S0) { + /* opposite direction */ + return taka_runge_kutta_reverse( + taka_runge_kutta_4a(map(taka_minus,base_replace(FF,[[X0,-X0]])),X0,Y,-S0,Ys,-T0,H)); + } + if (H >= T0-S0) { + A=taka_runge_kutta_4(FF,X0,Y,S0,Ys,T0-S0,0 | onestep=1); + }else{ + A=taka_runge_kutta_4(FF,X0,Y,S0,Ys,H,T0); + T=A[0]; + if (T0-T[0] > 0) { + B=taka_runge_kutta_4(FF,X0,Y,T[0],cdr(T),T0-T[0],0 | onestep=1); + T=B[0]; + A=cons(T,A); + } + } + return(A); +} +def runge_kutta_4(FF,X0,Y,S0,Ys,T0,H) { + return taka_runge_kutta_4a(FF,X0,Y,S0,Ys,T0,H); +} +def runge_kutta_4_linear(FF,X0,Y,S0,Ys,T0,H) { + return taka_runge_kutta_4a_linear(FF,X0,Y,S0,Ys,T0,H); +} + +def taka_runge_kutta_4a_linear(FF,X0,Y,S0,Ys,T0,H) { + if (T0 < S0) { + /* opposite direction */ + return taka_runge_kutta_reverse( + taka_runge_kutta_4a_linear(map(taka_minus,base_replace(FF,[[X0,-X0]])),X0,Y,-S0,Ys,-T0,H)); + } + if (H >= T0-S0) { + A=taka_runge_kutta_4_linear(FF,X0,Y,S0,Ys,T0-S0,0 | onestep=1); + }else{ + A=taka_runge_kutta_4_linear(FF,X0,Y,S0,Ys,H,T0); + T=A[0]; + if (T0-T[0] > 0) { + B=taka_runge_kutta_4_linear(FF,X0,Y,T[0],cdr(T),T0-T[0],0 | onestep=1); + T=B[0]; + A=cons(T,A); + } + } + return(A); +} + +/* equation of oscilation */ +def test4() { + A=runge_kutta_4([y1,-y0],x,[y0,y1],0,[1,0],3.14*2,0.1); + taka_plot_auto(A); + return(A); +} + +def test4b() { + A=taka_runge_kutta_4a([y1,-y0],x,[y0,y1],3.14,[-1,0],0,0.1); + taka_plot_auto(A); + return(A); +} +endmodule; + +import("taka_plot.rr")$ +pari(allocatemem,10^7)$ +module rktest; +localf re$ +localf re2$ +localf im$ +localf im2$ +localf tryA$ +localf tryA2$ +localf geq$ +localf test1$ + +def re(L) { + return map(re2,L); +} +def re2(P) { + return map(number_real_part,P); +} +def im(L) { + return map(im2,L); +} +def im2(P) { + return map(number_imaginary_part,P); +} + +def geq() { + L=x*(1-x)*dx^2+(c-(a+b+1)*x)*dx-a*b; + L=base_replace(L,[[a,1/2],[b,1/2],[c,1]]); + + L2 = -((c-(a+b+1)*x)*y2-a*b*y1)/(x*(1-x)); + L2=base_replace(L2,[[a,1/2],[b,1/2],[c,1]]); + return [ y2, L2]; +} +def tryA() { + LL = geq(); + A = taka_runge_kutta_4_adaptive( + LL, + x,[y1,y2], + 0.5+0.5*@i,[-6.78383-1.28991*@i, -1.51159-1.7935*@i], + (3-@i)*0.0005, 2.0); + taka_plot_auto(re(A)); + return A; +} +def tryA2() { + LL = geq(); + A = taka_runge_kutta_4( + LL, + x,[y1,y2], + 0.5+0.5*@i,[-6.78383-1.28991*@i, -1.51159-1.7935*@i], + (3-@i)*0.0005, 1.0); + taka_plot_auto(re(A)); + return A; +} + +def test1() { + A=tryA2(); + B=A[100]; + print(B); +/* cf. misc-2008/A2/misc/ellip2.m + p2 = 0.5+0.5*I --> [-6.78383-1.28991*@i, -1.51159-1.7935*@i]. + p2 = 0.8495+0.3835*I; +Print["-------------------"]; +Print[p2]; +Print[N[-2*Gamma[1/2]^2*Hypergeometric2F1[1/2,1/2,1,z] /. {z->p2}]] +Print[N[D[-2*Gamma[1/2]^2*Hypergeometric2F1[1/2,1/2,1,z],z] /. {z->p2}]] +*/ + print("math: [0.8495+0.3835*I, -7.64079 - 2.0799*I, -1.16364 - 4.14709*I] "); + print("It was Buggy tryA2() and tryA() --> fixed. see log of 1.8"); +} + +endmodule; + end$