share/contrib/rkf45/rkf45.mac

   1 /*
   2 --------------------------------------------------------------------------------
   3 Copyright (C) 2011 Panagiotis Papasotiriou
   4
   5 This software is released under the terms of the GNU General Public License.
   6 See <http://www.gnu.org/licenses/> for details.
   7 --------------------------------------------------------------------------------
   8      Version: 1
   9 Date created: September 12, 2011
  10  Last update: October 26, 2011
  11       Author: Panagiotis J. Papasotiriou
  12 --------------------------------------------------------------------------------
  13 Brief description:
  14
  15 rkf45 is a Maxima function for solving initial value problems with automatic
  16 step size and error control.
  17 This is an implementation of the Runge-Kutta-Fehlberg 4th-5th-order scheme.
  18 --------------------------------------------------------------------------------
  19 Syntax:
  20
  21 rkf45(ode,func,init,interval,options)
  22 rkf45([ode1,ode2,...],[func1,func2,...],[init1,init2,...],interval,options)
  23
  24 The first form solves a first-order differential equation, ode, with respect to
  25 the initial condition init, where func is the dependent variable and init is the
  26 value of the dependent variable at the initial point.
  27
  28 Similarly, the second form solves a system of first-order differential
  29 equations, ode1,ode2,..., with respect to the initial conditions
  30 init1,init2,..., where func1,func2,... are the dependent variables and
  31 init1,init2,... are the values of the dependent variables at the initial point.
  32
  33 Differential equation(s) should be given as expressions depending only on the
  34 independent and dependent variables, and should define the derivative of the
  35 dependent variable with respect to the independent variable. For instance, the
  36 differential equation y'(x)+(x+1)*y=0 should be given as -(x+1)*y.
  37 The argument "interval" should be a list of three elements. The first element
  38 identifies the independent variable, while the second and third elements are the
  39 initial and final values for the independent variable, for instance: [x,0,6].
  40 Initial value does not need to be less than final value, so an interval like
  41 [x,6,0] is also valid.
  42
  43 rkf45 accepts the following optional arguments:
  44
  45 *      full_solution: A Boolean. If set to true, a full list of the solution at
  46                       all intermediate points will be returned. If set to false,
  47                       only the solution at the last integration point is
  48                       returned. Default: true.
  49
  50 * absolute_tolerance: The desired absolute tolerance of the result. Default:
  51                       1e-6.
  52
  53 *     max_iterations: Maximum number of iterations. Default: 10000.
  54
  55 *            h_start: Initial integration step. Default: one 100th of the
  56                       integration interval, (interval[3]-interval[2])/100.
  57
  58 *             report: A Boolean. If set to true, rkf45 prints a report at
  59                       exit, giving details about the calculations done.
  60                       Default: false.
  61
  62 Integration step size ia selected automatically in such a way that the estimated
  63 local error is less than user-specified absolute tolerance.
  64 The result is returned as a list with n+1 columns, where n is the number of
  65 first-order differential equations. The first column contains the values of the
  66 independent variable selected by the algorithm. The second column contains the
  67 values of the first dependent variable at the corresponding value of the
  68 independent variable. Similarly, the third column contains the values of the
  69 second dependent variable at the corresponding value of the independent
  70 variable, and so on.
  71
  72 rkf45 can be used to solve moderately stiff initial value problems, although it
  73 is not designed for that purpose.
  74 --------------------------------------------------------------------------------
  75 Examples:
  76
  77 (1) A first-order differential equation, y'=x*(y-1)+3, with y(0)=-2:
  78     rkf45(x*(y-1)+3,y,-2,[x,0,4]) returns solution at selected points from x=0
  79     to x=4.
  80
  81 (2) A second-order differential equation, y''=x+y*y', with y(-1)=2, y'(-1)=0:
  82     rkf45([y2,x+y1*y2],[y1,y2],[2,0],[x,-1,4]) returns solution at selected
  83     points from x=-1 to x=4.
  84
  85 See rkf45.dem for more examples. Detailed documentation and several examples
  86 discussed in detail can be found at the accompanying document rkf45.pdf.
  87 --------------------------------------------------------------------------------
  88 */
  89 rkf45(odes,funcs,initial,interval,[options]):=block(
  90   [numer:true,atol,save_steps,maxit,show_report,xi,xc,yi,h,not_ok,
  91    k1,k2,k3,k4,k5,k6,trunc_error,min_trunc_error,h_optimal,estimated_errors,step_extrema,
  92    x_tiny,bad_steps:0,sol],
  93   /* Check interval for errors */
  94   if (not(listp(interval))) or length(interval)<3 then error("Error: interval should be a list of the independent variable, start and end times and optionally, a stepsize. Found",interval,"instead."),
  95   x_tiny:1e-7*interval[3],
  96   /* Set optional arguments */
  97   atol:assoc('absolute_tolerance,options,1e-6),
  98   save_steps:assoc('full_solution,options,true),
  99   maxit:assoc('max_iterations,options,10000),
 100   show_report:assoc('report,options,false),
 101   h:assoc('h_start,options,(interval[3]-interval[2])/100),
 102   min_trunc_error:assoc('min_trunc_error,options,1e-30),
 103   /* define setter of h_optimal */
 104   rkf45_h_optimal:lambda([h,atol,trunc_error],min(max(0.84089641525371*(atol/max(min_trunc_error,trunc_error))^0.25,0.1),4)*h),
 105   /* Convert arguments to lists, if necessary */
 106   if (not(listp(odes))) then odes:[odes],
 107   if (not(listp(funcs))) then funcs:[funcs],
 108   if (not(listp(initial))) then initial:[initial],
 109   /* coerce odes, interval and initial to floats */
 110   [odes,interval,initial]:float([odes,interval,initial]),
 111   /* Define right-hand-side function */
 112   local(f_rhs),
 113   block([vars:cons(interval[1],funcs)],
 114     define(funmake(f_rhs,vars),buildq([odes:odes,v:funmake(mode_declare,flatten(map(lambda([t],[t,'float]),vars)))],block(v,odes)))),
 115   translate(f_rhs),
 116   /* Initialize */
 117   xi:interval[2], yi:initial,
 118   /* Check initial values */
 119   if not every(map(numberp,initial)) then  error("Error: initial should be a list of numbers, but found",initial,"instead."),
 120   block([initial_values:apply(f_rhs,cons(xi,yi))],
 121     if not every(map(numberp,initial_values)) then error("Error: odes should evaluate to a list of numbers at initial, but found",initial_values,"instead.")),
 122   /* Set up the rest */
 123   estimated_errors:[1e30,-1e30],
 124   step_extrema:[1e30,-1e30],
 125   if save_steps then sol:[cons(xi,yi)],
 126   not_ok:true,
 127   /* Main loop */
 128   for i:1 thru maxit do (
 129     /* Compute solution at xi+h */
 130     xc:xi,
 131     k1:h*apply(f_rhs,cons(xi,yi)),
 132     k2:h*apply(f_rhs,cons(xi+0.25*h,yi+0.25*k1)),
 133     k3:h*apply(f_rhs,cons(xi+0.375*h,yi+0.09375*k1+0.28125*k2)),
 134     k4:h*apply(f_rhs,cons(xi+9.230769230769231e-1*h,yi+8.793809740555303e-1*k1
 135                                                       -3.277196176604461*k2
 136                                                       +3.320892125625854*k3)),
 137     k5:h*apply(f_rhs,cons(xi+h,yi+2.032407407407407*k1-8*k2+7.173489278752437*k3
 138                                  -2.058966861598441e-1*k4)),
 139     k6:h*apply(f_rhs,cons(xi+0.5*h,yi-2.962962962962963e-1*k1+2*k2
 140                                      -1.381676413255361*k3
 141                                      +4.529727095516569e-1*k4-0.275*k5)),
 142     /* Estimate local truncation error */
 143     trunc_error:lmax(abs(2.777777777777778e-3*k1-2.994152046783626e-2*k3
 144                         -2.919989367357788e-2*k4+0.02*k5+3.636363636363636e-2*k6
 145                         ))/abs(h),
 146     if (trunc_error<atol) then (
 147       /* Step accepted, do it */
 148       xi:xc+h,
 149       yi:yi+1.157407407407407e-1*k1+5.489278752436647e-1*k3
 150            +5.353313840155946e-1*k4-0.2*k5,
 151       if save_steps then sol:endcons(cons(xi,yi),sol),
 152       estimated_errors[1]:min(trunc_error,estimated_errors[1]),
 153       estimated_errors[2]:max(trunc_error,estimated_errors[2]),
 154       step_extrema[1]:min(h,step_extrema[1]),
 155       step_extrema[2]:max(h,step_extrema[2]),
 156       /* Return if final point is reached. */
 157       if abs(xi-interval[3])<=x_tiny then (
 158         not_ok:false,
 159         if not(save_steps) then sol:cons(xi,yi),
 160         if show_report then (
 161           print("------------------------------------------------------"),
 162           print("Info: rkf45:"),
 163           print("   Integration points selected:",i+1-bad_steps),
 164           print("    Total number of iterations:",i),
 165           print("           Bad steps corrected:",bad_steps),
 166           print("       Minimum estimated error:",estimated_errors[1]),
 167           print("       Maximum estimated error:",estimated_errors[2]),
 168           print("Minimum integration step taken:",step_extrema[1]),
 169           print("Maximum integration step taken:",step_extrema[2]),
 170           print("------------------------------------------------------")
 171         ),
 172         return()
 173       )
 174     ) else (
 175       /* Step is not accepted, must try again with optimal step */
 176       bad_steps:bad_steps+1
 177     ),
 178     /* Prepare next step */
 179     h_optimal:rkf45_h_optimal(h,atol,trunc_error),
 180     if h_optimal>0 then
 181       if xi+h_optimal<interval[3] then h:h_optimal else h:interval[3]-xi
 182     else if xi+h_optimal>interval[3] then h:h_optimal else h:interval[3]-xi
 183   ),
 184   /* Warn user if necessary */
 185   if not_ok then (
 186     print("Warning: rkf45: Integration stopped at x =",xi,"(stiff problem?)"),
 187     print("                Iterations limit has been reached. Check if differential"),
 188     print("                equations/initial conditions are given correctly, reduce"),
 189     print("                accuracy, and/or increase maximum number of steps.")
 190   ),
 191   /* Return solution */
 192   return(sol)
 193 )$
 194 /*----------------------------------------------------------------------------*/