path: root/src/algebra/numode.spad.pamphlet

\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}
<<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(try: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]
      i       : I

      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 (x + step.try - x2) > 0.0 then
            step.try := x2 - x
--find the correct scaling
         for i in 1..m repeat
            yscal(i) := abs(y(i)) + abs(step.try * 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.try := 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.try
      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
      i      : I
      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
     i    : I
     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
      i  : I
-- 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}
<<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.
@
<<*>>=
<<license>>

<<package NUMODE NumericalOrdinaryDifferentialEquations>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}