#### odeex.red3.1 KB Permalink History Raw

 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123 ``````% A SIMPLE PROGRAM FOR COMPUTING SOLUTIONS OF ODES BY TAYLOR SERIES. % Author: Andreas Strotmann . 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 <> 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 <> 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; ``````