Initial population.

1 files changed, 410 insertions, 0 deletions

diff --git a/src/algebra/numode.spad.pamphlet b/src/algebra/numode.spad.pamphlet
new file mode 100644
index 00000000..4de1cf27
--- /dev/null
+++ b/src/algebra/numode.spad.pamphlet
@@ -0,0 +1,410 @@
+\documentclass{article}
+\usepackage{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
+      iter    : I
+
+      eps  := 1.0/eps
+      for i in 1..m repeat
+         y(i)  := ystart(i)
+      for iter in 1..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
+--       outlist.0 := x::E
+--       outlist.1 := y(0)::E
+--       outlist.2 := y(1)::E
+--       output(blankSeparate(outlist)::E)
+--check to see if done
+         if (x-x2) >= 0.0 then
+            leave
+--next stepsize to use
+         step.try := step.next
+--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}