axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / NUMODE.spad

)abbrev package NUMODE NumericalOrdinaryDifferentialEquations
++ Author: Yurij Baransky
++ Date Created: October 90
++ Date Last Updated: October 90
++ Description:
++ This package is a suite of functions for the numerical integration of an
++ ordinary differential equation of n variables:\br
++ \tab{5}dy/dx = f(y,x)\tab{5}y is an n-vector\br
++ All the routines are based on a 4-th order Runge-Kutta kernel.
++ These routines generally have as arguments:\br
++ n, the number of dependent variables;\br
++ x1, the initial point;\br
++ h, the step size;\br
++ y, a vector of initial conditions of length n\br
++ which upon exit contains the solution at \spad{x1 + h};\br
++
++ \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.
++
++ In order of increasing complexity:\br
++ \tab{5}\spad{rk4(y,n,x1,h,derivs)} advances the solution vector to\br
++ \tab{5}\spad{x1 + h} and return the values in y.\br
++
++ \tab{5}\spad{rk4(y,n,x1,h,derivs,t1,t2,t3,t4)} is the same as\br
++ \tab{5}\spad{rk4(y,n,x1,h,derivs)} except that you must provide 4 scratch\br
++ \tab{5}arrays t1-t4 of size n.\br
++
++ \tab{5}Starting with y at x1, \spad{rk4f(y,n,x1,x2,ns,derivs)}\br
++ \tab{5}uses \spad{ns} fixed steps of a 4-th order Runge-Kutta\br
++ \tab{5}integrator to advance the solution vector to x2 and return\br
++ \tab{5}the values in y.  Argument x2, is the final point, and\br
++ \tab{5}\spad{ns}, the number of steps to take.
++
++ \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\br
++ \tab{5}\spad{error = MAX(ABS((y2steps(i) - y1step(i))/yscal(i)))}\br
++ and this is compared against \spad{eps}. If this is greater
++ than \spad{eps}, the step size is reduced accordingly to\br
++ \tab{5}\spad{hnew = 0.9 * hdid * (error/eps)**(-1/4)}\br
++ 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\br
++ \tab{5}\spad{hnext = 0.9 * hdid * (error/\spad{eps})**(-1/5)}\br
++ 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.
++
++ \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.
++
++ \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\br
++ \tab{5}\spad{yscal(i) = abs(y(i)) + abs(h*dydx(i)) + tiny}\br
++ 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.

NumericalOrdinaryDifferentialEquations() : SIG == CODE 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)

  SIG ==> with

    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
      ++ 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
      ++ 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
      ++ 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 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
      ++ 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}.

  CODE ==>  add

   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
         --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) )