\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra numode.spad} \author{Yurij Baransky} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package NUMODE NumericalOrdinaryDifferentialEquations} <>= )abbrev package NUMODE NumericalOrdinaryDifferentialEquations ++ Author: Yurij Baransky ++ Date Created: October 90 ++ Date Last Updated: October 90 ++ Basic Operations: ++ Related Constructors: ++ Also See: ++ AMS Classifications: ++ Keywords: ++ References: ++ Description: ++ This package is a suite of functions for the numerical integration of an ++ ordinary differential equation of n variables: ++ ++ \center{dy/dx = f(y,x)\space{5}y is an n-vector} ++ ++ \par All the routines are based on a 4-th order Runge-Kutta kernel. ++ These routines generally have as arguments: ++ n, the number of dependent variables; ++ x1, the initial point; ++ h, the step size; ++ y, a vector of initial conditions of length n which upon exit contains the solution at \spad{x1 + h}; ++ \spad{derivs}, a function which computes the right hand side of the ++ ordinary differential equation: \spad{derivs(dydx,y,x)} computes \spad{dydx}, ++ a vector which contains the derivative information. ++ ++ \par In order of increasing complexity:\begin{items} ++ ++ \item \spad{rk4(y,n,x1,h,derivs)} advances the solution vector to ++ \spad{x1 + h} and return the values in y. ++ ++ \item \spad{rk4(y,n,x1,h,derivs,t1,t2,t3,t4)} is the same as ++ \spad{rk4(y,n,x1,h,derivs)} except that you must provide 4 scratch ++ arrays t1-t4 of size n. ++ ++ \item Starting with y at x1, \spad{rk4f(y,n,x1,x2,ns,derivs)} ++ uses \spad{ns} fixed ++ steps of a 4-th order Runge-Kutta integrator to advance the ++ solution vector to x2 and return the values in y. ++ Argument x2, is the final point, and ++ \spad{ns}, the number of steps to take. ++ ++ \item \spad{rk4qc(y,n,x1,step,eps,yscal,derivs)} takes a 5-th order ++ Runge-Kutta step with monitoring ++ of local truncation to ensure accuracy and adjust stepsize. ++ The function takes two half steps and one full step and scales ++ the difference in solutions at the final point. If the error is ++ within \spad{eps}, the step is taken and the result is returned. ++ If the error is not within \spad{eps}, the stepsize if decreased ++ and the procedure is tried again until the desired accuracy is ++ reached. Upon input, an trial step size must be given and upon ++ return, an estimate of the next step size to use is returned as ++ well as the step size which produced the desired accuracy. ++ The scaled error is computed as ++ \center{\spad{error = MAX(ABS((y2steps(i) - y1step(i))/yscal(i)))}} ++ and this is compared against \spad{eps}. If this is greater ++ than \spad{eps}, the step size is reduced accordingly to ++ \center{\spad{hnew = 0.9 * hdid * (error/eps)**(-1/4)}} ++ If the error criterion is satisfied, then we check if the ++ step size was too fine and return a more efficient one. If ++ \spad{error > \spad{eps} * (6.0E-04)} then the next step size should be ++ \center{\spad{hnext = 0.9 * hdid * (error/\spad{eps})**(-1/5)}} ++ Otherwise \spad{hnext = 4.0 * hdid} is returned. ++ A more detailed discussion of this and related topics can be ++ found in the book "Numerical Recipies" by W.Press, B.P. Flannery, ++ S.A. Teukolsky, W.T. Vetterling published by Cambridge University Press. ++ Argument \spad{step} is a record of 3 floating point ++ numbers \spad{(try , did , next)}, ++ \spad{eps} is the required accuracy, ++ \spad{yscal} is the scaling vector for the difference in solutions. ++ On input, \spad{step.try} should be the guess at a step ++ size to achieve the accuracy. ++ On output, \spad{step.did} contains the step size which achieved the ++ accuracy and \spad{step.next} is the next step size to use. ++ ++ \item \spad{rk4qc(y,n,x1,step,eps,yscal,derivs,t1,t2,t3,t4,t5,t6,t7)} is the ++ same as \spad{rk4qc(y,n,x1,step,eps,yscal,derivs)} except that the user ++ must provide the 7 scratch arrays \spad{t1-t7} of size n. ++ ++ \item \spad{rk4a(y,n,x1,x2,eps,h,ns,derivs)} ++ is a driver program which uses \spad{rk4qc} to integrate n ordinary ++ differential equations starting at x1 to x2, keeping the local ++ truncation error to within \spad{eps} by changing the local step size. ++ The scaling vector is defined as ++ \center{\spad{yscal(i) = abs(y(i)) + abs(h*dydx(i)) + tiny}} ++ where \spad{y(i)} is the solution at location x, \spad{dydx} is the ++ ordinary differential equation's right hand side, h is the current ++ step size and \spad{tiny} is 10 times the ++ smallest positive number representable. ++ The user must supply an estimate for a trial step size and ++ the maximum number of calls to \spad{rk4qc} to use. ++ Argument x2 is the final point, ++ \spad{eps} is local truncation, ++ \spad{ns} is the maximum number of call to \spad{rk4qc} to use. ++ \end{items} NumericalOrdinaryDifferentialEquations(): Exports == Implementation where L ==> List V ==> Vector B ==> Boolean I ==> Integer E ==> OutputForm NF ==> Float NNI ==> NonNegativeInteger VOID ==> Void OFORM ==> OutputForm RK4STEP ==> Record(tryValue:NF, did:NF, next:NF) Exports ==> with --header definitions here rk4 : (V NF,I,NF,NF, (V NF,V NF,NF) -> VOID) -> VOID ++ rk4(y,n,x1,h,derivs) uses a 4-th order Runge-Kutta method ++ to numerically integrate the ordinary differential equation ++ {\em dy/dx = f(y,x)} of n variables, where y is an n-vector. ++ Argument y is a vector of initial conditions of length n which upon exit ++ contains the solution at \spad{x1 + h}, n is the number of dependent ++ variables, x1 is the initial point, h is the step size, and ++ \spad{derivs} is a function which computes the right hand side of the ++ ordinary differential equation. ++ For details, see \spadtype{NumericalOrdinaryDifferentialEquations}. rk4 : (V NF,I,NF,NF, (V NF,V NF,NF) -> VOID ,V NF,V NF,V NF,V NF) -> VOID ++ rk4(y,n,x1,h,derivs,t1,t2,t3,t4) is the same as ++ \spad{rk4(y,n,x1,h,derivs)} except that you must provide 4 scratch ++ arrays t1-t4 of size n. ++ For details, see \con{NumericalOrdinaryDifferentialEquations}. rk4a : (V NF,I,NF,NF,NF,NF,I,(V NF,V NF,NF) -> VOID ) -> VOID ++ rk4a(y,n,x1,x2,eps,h,ns,derivs) is a driver function for the ++ numerical integration of an ordinary differential equation ++ {\em dy/dx = f(y,x)} of n variables, where y is an n-vector ++ using a 4-th order Runge-Kutta method. ++ For details, see \con{NumericalOrdinaryDifferentialEquations}. rk4qc : (V NF,I,NF,RK4STEP,NF,V NF,(V NF,V NF,NF) -> VOID) -> VOID ++ rk4qc(y,n,x1,step,eps,yscal,derivs) is a subfunction for the ++ numerical integration of an ordinary differential equation ++ {\em dy/dx = f(y,x)} of n variables, where y is an n-vector ++ using a 4-th order Runge-Kutta method. ++ This function takes a 5-th order Runge-Kutta step with monitoring ++ of local truncation to ensure accuracy and adjust stepsize. ++ For details, see \con{NumericalOrdinaryDifferentialEquations}. rk4qc : (V NF,I,NF,RK4STEP,NF,V NF,(V NF,V NF,NF) -> VOID ,V NF,V NF,V NF,V NF,V NF,V NF,V NF) -> VOID ++ rk4qc(y,n,x1,step,eps,yscal,derivs,t1,t2,t3,t4,t5,t6,t7) is a subfunction for the ++ numerical integration of an ordinary differential equation ++ {\em dy/dx = f(y,x)} of n variables, where y is an n-vector ++ using a 4-th order Runge-Kutta method. ++ This function takes a 5-th order Runge-Kutta step with monitoring ++ of local truncation to ensure accuracy and adjust stepsize. ++ For details, see \con{NumericalOrdinaryDifferentialEquations}. rk4f : (V NF,I,NF,NF,I,(V NF,V NF,NF) -> VOID ) -> VOID ++ rk4f(y,n,x1,x2,ns,derivs) uses a 4-th order Runge-Kutta method ++ to numerically integrate the ordinary differential equation ++ {\em dy/dx = f(y,x)} of n variables, where y is an n-vector. ++ Starting with y at x1, this function uses \spad{ns} fixed ++ steps of a 4-th order Runge-Kutta integrator to advance the ++ solution vector to x2 and return the values in y. ++ For details, see \con{NumericalOrdinaryDifferentialEquations}. Implementation ==> add --some local function definitions here rk4qclocal : (V NF,V NF,I,NF,RK4STEP,NF,V NF,(V NF,V NF,NF) -> VOID ,V NF,V NF,V NF,V NF,V NF,V NF) -> VOID rk4local : (V NF,V NF,I,NF,NF,V NF,(V NF,V NF,NF) -> VOID ,V NF,V NF,V NF) -> VOID import OutputPackage ------------------------------------------------------------ rk4a(ystart,nvar,x1,x2,eps,htry,nstep,derivs) == y : V NF := new(nvar::NNI,0.0) yscal : V NF := new(nvar::NNI,1.0) dydx : V NF := new(nvar::NNI,0.0) t1 : V NF := new(nvar::NNI,0.0) t2 : V NF := new(nvar::NNI,0.0) t3 : V NF := new(nvar::NNI,0.0) t4 : V NF := new(nvar::NNI,0.0) t5 : V NF := new(nvar::NNI,0.0) t6 : V NF := new(nvar::NNI,0.0) step : RK4STEP := [htry,0.0,0.0] x : NF := x1 tiny : NF := 10.0**(-(digits()+1)::I) m : I := nvar outlist : L OFORM := [x::E,x::E,x::E] eps := 1.0/eps for i in 1..m repeat y(i) := ystart(i) iter : I := 1 while iter <= nstep repeat --compute the derivative derivs(dydx,y,x) --if overshoot, the set h accordingly if positive?(x + step.tryValue - x2) then step.tryValue := x2 - x --find the correct scaling for i in 1..m repeat yscal(i) := abs(y(i)) + abs(step.tryValue * dydx(i)) + tiny --take a quality controlled runge-kutta step rk4qclocal(y,dydx,nvar,x,step,eps,yscal,derivs ,t1,t2,t3,t4,t5,t6) x := x + step.did --check to see if done if (x-x2) >= 0.0 then leave --next stepsize to use step.tryValue := step.next iter := iter + 1 --end nstep repeat if iter = (nstep+1) then output("ode: ERROR ") outlist.1 := nstep::E outlist.2 := " steps to small, last h = "::E outlist.3 := step.did::E output(blankSeparate(outlist)) output(" y= ",y::E) for i in 1..m repeat ystart(i) := y(i) ---------------------------------------------------------------- rk4qc(y,n,x,step,eps,yscal,derivs) == t1 : V NF := new(n::NNI,0.0) t2 : V NF := new(n::NNI,0.0) t3 : V NF := new(n::NNI,0.0) t4 : V NF := new(n::NNI,0.0) t5 : V NF := new(n::NNI,0.0) t6 : V NF := new(n::NNI,0.0) t7 : V NF := new(n::NNI,0.0) derivs(t7,y,x) eps := 1.0/eps rk4qclocal(y,t7,n,x,step,eps,yscal,derivs,t1,t2,t3,t4,t5,t6) -------------------------------------------------------- rk4qc(y,n,x,step,eps,yscal,derivs,t1,t2,t3,t4,t5,t6,dydx) == derivs(dydx,y,x) eps := 1.0/eps rk4qclocal(y,dydx,n,x,step,eps,yscal,derivs,t1,t2,t3,t4,t5,t6) -------------------------------------------------------- rk4qclocal(y,dydx,n,x,step,eps,yscal,derivs ,t1,t2,t3,ysav,dysav,ytemp) == xsav : NF := x h : NF := step.tryValue fcor : NF := 1.0/15.0 safety : NF := 0.9 grow : NF := -0.20 shrink : NF := -0.25 errcon : NF := 0.6E-04 --(this is 4/safety)**(1/grow) hh : NF errmax : NF m : I := n -- for i in 1..m repeat dysav(i) := dydx(i) ysav(i) := y(i) --cut down step size till error criterion is met repeat --take two little steps to get to x + h hh := 0.5 * h rk4local(ysav,dysav,n,xsav,hh,ytemp,derivs,t1,t2,t3) x := xsav + hh derivs(dydx,ytemp,x) rk4local(ytemp,dydx,n,x,hh,y,derivs,t1,t2,t3) x := xsav + h --take one big step get to x + h rk4local(ysav,dysav,n,xsav,h,ytemp,derivs,t1,t2,t3) --compute the maximum scaled difference errmax := 0.0 for i in 1..m repeat ytemp(i) := y(i) - ytemp(i) errmax := max(errmax,abs(ytemp(i)/yscal(i))) --scale relative to required accuracy errmax := errmax * eps --update integration stepsize if (errmax > 1.0) then h := safety * h * (errmax ** shrink) else step.did := h if errmax > errcon then step.next := safety * h * (errmax ** grow) else step.next := 4 * h leave --make fifth order with 4-th order error estimate for i in 1..m repeat y(i) := y(i) + ytemp(i) * fcor -------------------------------------------- rk4f(y,nvar,x1,x2,nstep,derivs) == yt : V NF := new(nvar::NNI,0.0) dyt : V NF := new(nvar::NNI,0.0) dym : V NF := new(nvar::NNI,0.0) dydx : V NF := new(nvar::NNI,0.0) ynew : V NF := new(nvar::NNI,0.0) h : NF := (x2-x1) / (nstep::NF) x : NF := x1 j : I -- start integrating for i in 1..nstep repeat derivs(dydx,y,x) rk4local(y,dydx,nvar,x,h,y,derivs,yt,dyt,dym) x := x + h -------------------------------------------------------- rk4(y,n,x,h,derivs) == t1 : V NF := new(n::NNI,0.0) t2 : V NF := new(n::NNI,0.0) t3 : V NF := new(n::NNI,0.0) t4 : V NF := new(n::NNI,0.0) derivs(t1,y,x) rk4local(y,t1,n,x,h,y,derivs,t2,t3,t4) ------------------------------------------------------------ rk4(y,n,x,h,derivs,t1,t2,t3,t4) == derivs(t1,y,x) rk4local(y,t1,n,x,h,y,derivs,t2,t3,t4) ------------------------------------------------------------ rk4local(y,dydx,n,x,h,yout,derivs,yt,dyt,dym) == hh : NF := h*0.5 h6 : NF := h/6.0 xh : NF := x+hh m : I := n -- first step for i in 1..m repeat yt(i) := y(i) + hh*dydx(i) -- second step derivs(dyt,yt,xh) for i in 1..m repeat yt(i) := y(i) + hh*dyt(i) -- third step derivs(dym,yt,xh) for i in 1..m repeat yt(i) := y(i) + h*dym(i) dym(i) := dyt(i) + dym(i) -- fourth step derivs(dyt,yt,x+h) for i in 1..m repeat yout(i) := y(i) + h6*( dydx(i) + 2.0*dym(i) + dyt(i) ) @ \section{License} <>= --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- - Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- -- - Redistributions in binary form must reproduce the above copyright -- notice, this list of conditions and the following disclaimer in -- the documentation and/or other materials provided with the -- distribution. -- -- - Neither the name of The Numerical ALgorithms Group Ltd. nor the -- names of its contributors may be used to endorse or promote products -- derived from this software without specific prior written permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS --IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED --TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A --PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER --OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, --EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, --PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @ <<*>>= <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}