Return to taka_runge_kutta.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

taka_runge_kutta.rr : the fourth order adaptive Runge-Kutta method.
taka_pfp.rr         : Evaluating p F p-1 by using the fourth order adaptive
                      Runge-Kutta method. (work in progress)
                      As for algorithmic aspects, see Numerical Recipes.

/* $Id$ */
/* 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 = 0$ /* Default 1,  H = 2*H */

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*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(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*0.5;
    }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);
  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*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(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);
    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;
      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)
          && Taka_Runge_kutta_Make_Larger) {
        H = 2*H;
        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() {
   /* 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);
   taka_plot_auto(T);
   print([T[0][0],eval(exp(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$