reduce-historical/r38/packages/numeric/rungeku.red

module rungeku;  % Numeric solution of ODE's with Runge-Kutta.
   % Version 2: supporting ODE systems

% Author: H. Melenk, ZIB, Berlin

% Copyright (c): ZIB Berlin 1993, all rights resrved

fluid '(!*noequiv accuracy!*);
global '(iterations!* !*trnumeric);

put ('num_odesolve,'psopfn,'rungekuttaeval);

symbolic procedure rungekuttaeval u;
     % interface function;
  begin scalar e,f,x,y,sx,sy,en,d,v,w,q;
    u := for each x in u collect reval x;
    u := accuracycontrol(u,20,6);

     % equations
    e :=car u; u :=cdr u;
    e := if eqcar(e,'list) then cdr e else {e};

     % dependent variable
    q :=car u; u :=cdr u;
    q := if eqcar(q,'list) then cdr q else {q};
    for each yy in q do
    <<if not eqcar(yy,'equal) or not idp cadr yy
     then typerr(yy,"expression `dep. variable=starting value'");
     sy:=caddr yy.sy; y := cadr yy.y;
    >>;
    sy:=reversip sy; y:=reversip y;
    if length sy neq length e then
      rederr "number of equations and variables don't correspond";

      % independent variable
    x :=car u; u :=cdr u;
    if not eqcar(x,'equal) or not idp cadr x
       or null (w:=revalnuminterval(caddr x,t))
        then typerr(x,"expression `indep. variable=interval'");
    sx:=car w; en:=cadr w; x := cadr x;

      % convert expressions to explicit ODE system.
    q := for each yy in y collect {'df,yy,x};
    v:=cdr solveeval list('list. e, 'list . q);
    if cdr v then
      ruksyserr(nil,"cannot convert to explicit ODE system");
    f := for each d in q collect rukufindrhs(d,v);
    return rungekutta1(f,x,y,sx,en,sy);
end;

symbolic procedure rukufindrhs(d,e);
  % Find the righthand side for d in list e.
  if atom e then
      ruksyserr(d," cannot be moved to lhs of system") else
  if eqcar(car e,'equal) and cadr car e = d then
           reval caddr car e else
  if pairp e and eqcar(car e,'list) then  rukufindrhs(d,cdar e)
  else rukufindrhs(d,cdr e);

symbolic procedure ruksyserr(u,m);
  <<if u then writepri(mkquote u,'first);
    writepri(mkquote m,'last);
    rederr "Runge-Kutta failed";
  >>;

symbolic procedure rungekutta1(f,x,y,sx,ex,sy);
    %preparations for rungekutta iteration.
 (begin scalar acc,r,oldmode,!*noequiv;
    integer prec;
    if not memq(dmode!*,'(!:rd!: !:cr))then
       <<oldmode:=t; setdmode('rounded,t)>>;
    !*noequiv:=t; prec:=precision 0;
    sx := force!-to!-dm numr simp sx;
    ex := force!-to!-dm numr simp ex;
    sy := for each z in sy collect force!-to!-dm numr simp z;
    acc := !:!:quotient(1,expt(10,accuracy!*));
    if !*trnumeric then prin2t "starting rungekutta iteration";
    r := rungekutta2(f,x,y,sx,ex,sy,acc);
    if oldmode then setdmode('rounded,nil);
    precision prec;
    if null r then rederr "no solution found";
    return 'list.r;
  end) where !*roundbf=!*roundbf;

symbolic procedure rungekutta2(f,xx,yy,xs,xe,ys,acc);
  % Algorithm for numeric ODE solution
  % f(xx,yy): rhs;
  % x:     independent variable;                    
  % y:     dependent variable; 
  % s:     starting point;
  % e:     endpoint;                                
  % acc:   required accuracy
  % f,yy,ys are vectors (lists).
 dm!:
  begin scalar h,h2,h4,x,y,r,w1,w2,vars,dir;
    integer count,st;
    vars := xx.yy;
    x:=xs; y:=ys;
    h:=(xe - xs)/iterations!*;
     dir := h>0; st:= iterations!*;
     % compute initial stepsize
  adapt:
    h2:=h/2; h4:=h/4;
     % full stepsize.
    w1:=rungekuttastep(f,vars,x,y,h,h2);
     % same with half stepsize.
    w2:=rungekuttastep(f,vars,x,y,h2,h4);
    w2:=rungekuttastep(f,vars,x+h2,w2,h2,h4);
      % abs(w2 - w1)
    if normlist list!-list(w2,w1) > acc then
    <<h:=h2; st:=st+st; goto adapt>>;
    if !*trnumeric and not(st=iterations!*) then
      <<prin2
         "*** RungeKutta increased number of steps to ";
        prin2t st>>;
  loop:
    if (dir and x>xe) or (not dir and x<xe) then goto finis;
    r:=('list.x.y) . r;
    count:=add1 count;
    y:= rungekuttastep(f,vars,x,y,h,h2);
      % Advance solution.
    x:=x+h;
    goto loop;
  finis:
    r:= ('list.xx.yy).('list.xs.ys). rungekuttares(reversip r,st);
    return r;
  end;

symbolic procedure rungekuttares(l,st);
   % eliminate intermediate points.
     if st=iterations!* then l else
     << for i:=1:iterations!* collect
         <<for j:=1:m do l:= cdr l; car l>>
     >> where m=st/iterations!*;

symbolic procedure rungekuttastep(f,vars,x,y,h,h2);
  dm!:
   begin scalar d1,d2,d3,d4;
             %  f(x,y)
    d1:= list!-evaluate(f,vars,x.y);
             %  f(x+h2,y+h2*d1) 
    d2:= list!-evaluate(f,vars,
          (x+h2). list!+list(y,scal!*list(h2,d1)));
             %  f(x+h2,y+h2*d2)
    d3:= list!-evaluate(f,vars,
          (x+h2). list!+list(y,scal!*list(h2,d2)));
             %  f(x+h,y+h*d3)
    d4:= list!-evaluate(f,vars,
          (x+h).list!+list(y,scal!*list(h,d3)));
      % find y(x+h) = y+h*(d1 + 2*d2 + 2*d3 + d4)/6.
    y:=list!+list(y,
       scal!*list(!:!:quotient(h,6),
         list!+list(d1,
          list!+list(scal!*list(2,d2),
           list!+list(scal!*list(2,d3),
            d4)))));
    return y;
  end;

endmodule;

end;