/* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_runge_kutta.rr,v 1.4 2003/03/11 02:00:53 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$
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]);
}
}
}
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;
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);
K3 = newvect(N);
K4 = 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+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_save_data) {
Ans = cons(cons(Xk,vtol(Yk)),Ans);
}
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;
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 (H > 0) {
if (Xk > X1) break;
}else{
if (Xk < X1) 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);
}
Loaded_taka_runge_kutta=1$
end$