/* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_runge_kutta.rr,v 1.25 2017/07/09 23:57:36 takayama Exp $ */
/* From misc/200205/runge-kutta.rr */
/* They have not yet been registered in names.rr */
#define DEVAL(a) eval(a)
Taka_Runge_kutta_adapted0 = 0$
Taka_Runge_kutta_epsilon = 0.1$
Taka_Runge_kutta_H_Upper_Bound = 0.2$
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$
extern Tk_rk_deep_eval$
Tk_rk_deep_eval=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;
Beta = 0.5;
P = 1; Q = 1;
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 (type(F) == 5) {
N = size(F)[0];
}else{
N = length(F);
}
if (type(Y0) != 5) {
Y0 = newvect(N,Y0);
}
Yk = Y0;
K1 = newvect(N);
K2 = newvect(N);
Yk1 = newvect(N);
Xk = X0;
while (Xk < X1) {
taka_runge_kutta_replace(K1,F,Y,N,X,Xk,Yk);
taka_runge_kutta_replace(K2,F,Y,N,X,Xk+P*H,Yk+Q*H*K1);
Yk1 = Yk+H*Alpha*K1+H*Beta*K2;
if (Taka_Runge_kutta_save_data) {
Ans = cons(cons(Xk,vtol(Yk)),Ans);
}
print([Xk,Yk[0]]);
if (Taka_Runge_kutta_graphic0) glib_line(Xk,Yk[0],Xk+H,Yk1[0]);
Xk += H;
Yk = Yk1;
}
return Ans;
}
def taka_runge_kutta_2_test() {
/* Equation of oscilations */
F = newvect(2,[y2,-y1]);
Y = [y1,y2];
taka_runge_kutta_2(F,x,Y,0,[1,0],0.1,15);
}
def taka_runge_kutta_replace(V,F,Y,N,X,Xk,Rule_vector) {
for (I=0; I<N; I++) {
V[I] = subst(F[I],X,Xk);
for (J=0; J<N; J++) {
V[I] = subst(V[I],Y[J],Rule_vector[J]);
}
V[I] = eval(V[I]*exp(0));
}
}
def taka_runge_kutta_abs(V) {
if (type(V) != 5 && type(V) != 4) { /* not a vector */
if (ntype(V) == 4) { /* complex number */
return V*conj(V);
}else{
return(V*V);
}
}
if (type(V) == 5) N = size(V)[0];
if (type(V) == 4) N = length(V);
S = 0;
for (I=0; I<N; I++) {
if (ntype(V[I]) == 4) /* complex number */
S += V[I]*conj(V[I]);
else
S += V[I]*V[I];
}
return S;
}
def taka_runge_kutta_4(F,X,Y,X0,Y0,H,X1) {
/* N : rank of the ODE. */
extern Taka_Runge_kutta_adapted0, Taka_Runge_kutta_epsilon,
Taka_Runge_kutta_graphic0, Taka_Runge_kutta_yrange,
Taka_Runge_kutta_save_data, Taka_Runge_kutta_debug;
OneStep = getopt(onestep);
if (type(OneStep) <= 0) OneStep = 0; else OneStep = 1;
if (OneStep) X1=X0+2*H;
if ((H<0) && (X1-X0)>0) 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];
}else{
N = length(F);
}
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(K1,F,Y,N,X,Xk,Yk);
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);
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_adaptive(F,X,Y,X0,Y0,H,X1) {
/* N : rank of the ODE. */
extern Taka_Runge_kutta_epsilon,
Taka_Runge_kutta_graphic0, Taka_Runge_kutta_yrange,
Taka_Runge_kutta_save_data,
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();
glib_window(X0,Y0[0]-Taka_Runge_kutta_yrange,X1,Y0[0]+Taka_Runge_kutta_yrange);
}
if (type(F) == 5) {
N = size(F)[0];
}else{
N = length(F);
}
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);
Yk2 = newvect(N);
Yk3 = newvect(N);
Xk = X0;
while (true) {
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*(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*(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);
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 ((H<Taka_Runge_kutta_H_Upper_Bound) && Taka_Runge_kutta_Make_Larger) {
H = (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]]);
if (Taka_Runge_kutta_save_data) {
Ans = cons(cons(Xk,vtol(Yk)),Ans);
}
if (Taka_Runge_kutta_graphic0) glib_line(Xk,Yk[0],Xk+H,Yk1[0]);
}
}
return Ans;
}
/* load("glib"); load("taka_plot.rr"); to execute the functions below. */
def taka_runge_kutta_4_a_test() {
/* exponential function */
F = newvect(1,[y1]);
Y = [y1];
T = taka_runge_kutta_4_adaptive(F,x,Y,0,[1],0.1,5);
taka_plot_auto(T);
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]);
Y = [y1,y2];
T=taka_runge_kutta_4(F,x,Y,0,[1,0],0.1,15);
print(T);
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,Inhom) {
extern Tk_rk_deep_eval;
if (Tk_rk_deep_eval) {
V1=base_replace(F,[[X,Xk]]);
}else{
V1=base_replace_n(F,[[X,Xk]]);
}
V1=V1*Rule_vector;
for (I=0; I<N; I++) {
if (Tk_rk_deep_eval)
V[I] = eval(V1[I]*exp(0));
else
V[I] = V1[I];
}
if (Inhom != 0) {
/* printf("X=%a, Xk=%a\n",X,Xk); */
if (Tk_rk_deep_eval) {
V1=map(eval,base_replace(Inhom,[[X,Xk]]));
}else{
V1=base_replace_n(Inhom,[[X,Xk]]);
}
for (I=0; I<N; I++) V[I] += V1[I];
}
}
/* Y is a dummy */
def taka_runge_kutta_4_linear(F,X,Y,X0,Y0,H,X1) {
/* N : rank of the ODE. */
extern Taka_Runge_kutta_adapted0, Taka_Runge_kutta_epsilon,
Taka_Runge_kutta_graphic0, Taka_Runge_kutta_yrange,
Taka_Runge_kutta_save_data, Taka_Runge_kutta_debug,
Tk_rk_deep_eval;
if (type(getopt(number_eval))>0) {
Neval=1; Tk_rk_deep_eval=1;
} else Neval=0;
//printf("Neval=%a\n",Neval);
OneStep = getopt(onestep);
if (type(OneStep) <= 0) OneStep = 0; else OneStep = 1;
if (OneStep) X1=X0+2*H;
Inhom = getopt(inhom);
if (type(Inhom) < 0) Inhom = 0;
/* else printf("inhom=%a, onestep=%a\n",Inhom,OneStep); */ /* for debug */
if (type(Inhom) == 4) {
Inhom = newvect(length(Inhom),Inhom);
}
if ((H<0) && (X1-X0)>0) 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,Inhom);
taka_runge_kutta_replace_linear(K2,F,Y,N,X,Xk+H*(1/2),Yk+K1*(1/2)*H,Inhom);
taka_runge_kutta_replace_linear(K3,F,Y,N,X,Xk+H*(1/2),Yk+K2*(1/2)*H,Inhom);
taka_runge_kutta_replace_linear(K4,F,Y,N,X,Xk+H,Yk+K3*H,Inhom);
Yk1 = Yk+H*(K1/6+K2/3+K3/3+K4/6);
if (Neval) {Yk=number_eval(Yk); Yk1 = number_eval(Yk1);}
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) {
Opt = getopt();
return taka_runge_kutta_4a_linear(FF,X0,Y,S0,Ys,T0,H | option_list=Opt);
}
def taka_runge_kutta_4a_linear(FF,X0,Y,S0,Ys,T0,H) {
Opt = getopt();
if (T0 < S0) {
/* opposite direction */
if (H<0) H=-H;
Opt = [];
if (type(getopt(number_eval)) >= 0) Opt=[["number_eval",getopt(number_eval)]];
if (type(getopt(onestep)) >= 0) Opt=[["onestep",getopt(onestep)]];
if (type(getopt(inhom)) >= 0) Inhom = getopt(inhom); else Inhom = 0;
if (Inhom != 0) {
Inhom = map(taka_minus,base_replace(Inhom,[[X0,-X0]]));
Opt = append(Opt,[["inhom",Inhom]]);
}
return taka_runge_kutta_reverse(
taka_runge_kutta_4a_linear(map(taka_minus,base_replace(FF,[[X0,-X0]])),X0,Y,-S0,Ys,-T0,H | option_list=Opt));
}
if (H >= T0-S0) {
A=taka_runge_kutta_4_linear(FF,X0,Y,S0,Ys,T0-S0,0 | option_list=append(Opt,[["onestep",1]]));
}else{
A=taka_runge_kutta_4_linear(FF,X0,Y,S0,Ys,H,T0 | option_list=Opt);
T=A[0];
if (T0-T[0] > 0) {
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]]));
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$
localf rcheck1$
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");
}
def rcheck1() {
/* f'=f/z+z*cos(z) where f=z*sin(z) */
Tk_rk_deep_eval=1;
Step = 0.01;
X0 = 3.14/2;
X = 3.14;
Iv=[X0*1];
Eq=[[1/z]]; Inhom = [z*cos(z)];
A=tk_rk.runge_kutta_4_linear(Eq,z,[],X0,Iv,X,Step | inhom=Inhom, onestep=1);
print(A);
A=tk_rk.runge_kutta_4_linear(Eq,z,[],X0,Iv,X,Step | inhom=Inhom);
return A;
}
endmodule;
end$