File r35/lib/odeex.red artifact c3908f286d part of check-in b5833487d7

%  A SIMPLE PROGRAM FOR COMPUTING SOLUTIONS OF ODES BY TAYLOR SERIES.

% Author: Andreas Strotmann <strotmann@rrz.uni-koeln.de>.

algebraic;

% 1. The simplest case.

% Compute the first N terms of the Taylor series of the solution of
% the explicit ordinary first order differential equation
% y' = f(x,y)
% in a neighborhood of x0 *if* f is holomorphic in x and y at (x0,y0).

PROCEDURE detaylor(f,x,y,x0,y0,N);

begin scalar wj, pot;
  wj:=f; pot:=x-x0;
  return ( y0+sub({x=x0,y=y0},f)*(x-x0) +
            for j:=2:n sum
              << wj:= df(wj,x)+f*df(wj,y);
                 pot := pot*(x-x0)/j;
               sub({x=x0,y=y0},wj)*pot>>);
end;

% Example: y'=xy

detaylor(x*y,x,y,0,1,5);

% 2. The general case.

% Vectors (= systems of ODEs) are encoded as lists.

% 2.1 Auxiliaries.

infix lplusl;

precedence lplusl,+; 

procedure x lplusl y; % vector + vector
  begin scalar auxy;
    auxy:= y;
    return
     foreach xi in x collect <<auxy:= rest auxy; s>>
        where s= first auxy+ xi;
  end;

infix ltimesl;

precedence ltimesl,*;

procedure x ltimesl y; % vector * vector -> scalar
  begin scalar auxy;
    auxy:= y;
    return
     foreach xi in x sum <<auxy:= rest auxy; s>> where s=first auxy* xi;
  end;

infix ltimess;

precedence ltimess,*;

procedure x ltimess y; % vector * scalar  -> vector
  foreach xi in x collect y*xi;

% 2.2 The central procedure.

% Compute the first N terms of the Taylor series of the solution of
% the initial value problem
% (y1,...,yn)'=(f1(x,y1,...,yn), ... , fn(x,y1,...,yn))
%   such that  y1(x0)=y10, ..., yn(x0)=yn0
% for a system of explicit ordinary first order differential equations
% in a neighborhood of x0 *if* f is holomorphic in x and all the yi at
% (x0, y10,....,yn0).
%
% Input format:  flis={f1,...,fn},
%                Anfangswerte={x=x0, y1=y10,..., yn=yn0}
%
% NOTE: none of the yi may DEPEND on x (i.e., be symbols declared to
%       do so).
%       The yi MUST be symbols so DF can handle them.

procedure odetaylor(flis,Anfangswerte,N);
begin scalar pot,x,y,x0,y0,wj,res;
 % Split args (see comment above for format):
 x:=  lhs first Anfangswerte;
 x0:= rhs first Anfangswerte;
 y:=  for each gl in rest Anfangswerte collect lhs gl;
 y0:= for each gl in rest Anfangswerte collect rhs gl;
% Initialisations (= degree one of the taylor polynomial)
res:= y0 lplusl (sub(Anfangswerte,flis) ltimess (x-x0));
pot:= x-x0;
 wj:= flis;
% Main loop:
for j:=2:n do
   << wj:= foreach wij in wj
             collect df(wij,x) + (flis ltimesl
                                    foreach yk in y  % one row of the
                                                     % Jacobian
                                      collect df(wij,yk));  %of wj wrt y
% The above DFs should be PARTDFs, really. In REDUCE 3.4, maybe they
%    can...
     pot := pot*(x-x0)/j;
     res:= res lplusl (sub(Anfangswerte,wj) ltimess pot);
           %should be sub...
   >>;
 % DONE:
 return res;
end;

% Examples:

factor x;

% y''=-y.

odetaylor({yprime,-y}, {x=0,y=0,yprime=1}, 4);

% And something wild just for fun:

odetaylor({sin y2, cos(x*y1*y2)}, {x=0,y1=pi/2, y2=pi*7/4}, 4);

end;