Skip to content
Browse files

Commit package rkf45 (Runge-Kutta-Fehlberg 4th-5th order numerical o.…

…d.e. solver),

by Panagiotis Papasotiriou. Retrieved from: https://sites.google.com/site/pjpapasot/maxima/libraries/rkf45
  • Loading branch information...
1 parent a7a1015 commit c2c9482b316973c2bd636962d32ad6c439b7678e robert_dodier committed Apr 1, 2012
View
68 share/contrib/rkf45/rkf45.dem
@@ -0,0 +1,68 @@
+/* ---------------------------------------------------------------------------*/
+/* Demo file for rkf45. */
+/* ---------------------------------------------------------------------------*/
+load('rkf45)$
+/* ---------------------------------------------------------------------------*/
+( disp("","One differential equation of first order."),
+ sol:rkf45(-3*x*y^2+1/(x^3+1),y,0,[x,0,5],report=true),
+ plot2d([discrete,sol],[style,[lines,2]])
+)$
+/* ---------------------------------------------------------------------------*/
+( disp("","An initial value problem with threshold effect."),
+ sol:makelist(rkf45(equ,f,0,[t,0,100]),equ,
+ makelist(s-1.51*f+3.03*f^2/(1+f^2),s,[0.206,0.204,0.202])),
+ plot2d(makelist([discrete,s],s,sol),[style,[lines,2]],[xlabel,"t"],
+ [ylabel,"f"],[legend,"s=0.206","s=0.204","s=0.202"],
+ [gnuplot_preamble,"set key left"])
+)$
+/* ---------------------------------------------------------------------------*/
+( disp("","A system of two first-order differential equations."),
+ k1:0.4*8.8/62/0.03, k2:0.5*8.8/139/0.2/0.003, k3:0.4*8.8/139/0.2/0.003,
+ sol:rkf45([k1*(C-L),k2*(32-C)+k3*(L-C)],[L,C],[150,150],[t,0,2],report=true),
+ plot2d([[discrete,map(lambda([u],part(u,[1,2])),sol)],
+ [discrete,map(lambda([u],part(u,[1,3])),sol)]],
+ [style,[lines,2]],[xlabel,"time (hours)"],[ylabel,"temperature (F)"],
+ [legend,"liquid, L(t)","container, C(t)"])
+)$
+/* ---------------------------------------------------------------------------*/
+( disp("","A second-order differential equation: the van der Pol equation."),
+ equ:'diff(x,t,2)+4*(x^2-1)*'diff(x,t)+x=0,
+ equ2:['diff(x1,t)=x2,ev(solve(equ,'diff(x,t,2))[1],'diff(x,t,2)='diff(x2,t),
+ 'diff(x,t)=x2,x=x1)],
+ equ2:map(rhs,equ2),
+ sol:rkf45(equ2,[x1,x2],[0.75,0],[t,0,20],report=true),
+ plot2d([[discrete,map(lambda([u],part(u,[1,2])),sol)],
+ [discrete,map(lambda([u],part(u,[1,3])),sol)]],
+ [style,[lines,2]],[xlabel,"t"],[ylabel,"x, x'"],
+ [legend,"x(t)","x'(t)"])
+)$
+/* ---------------------------------------------------------------------------*/
+( disp("","A system of two second-order differential equations:",
+ "the double pendulum."),
+ m1:1, m2:1.5, l1:0.4, l2:0.6, g:9.81,
+ d2th1dt2:(-g*(2*m1+m2)*sin(th1)-m2*g*sin(th1-2*th2)
+ -2*sin(th1-th2)*m2*('diff(th2,t)^2*l2
+ +'diff(th1,t)^2*l1*cos(th1-th2)))
+ /(l1*(2*m1+m2-m2*cos(2*th1-2*th2))),
+ d2th2dt2:(2*sin(th1-th2)*('diff(th1,t)^2*l1*(m1+m2)+g*(m1+m2)*cos(th1)
+ +'diff(th2,t)^2*l2*m2*cos(th1-th2)))
+ /(l2*(2*m1+m2-m2*cos(2*th1-2*th2))),
+ equs:ev([omega1,omega2,d2th1dt2,d2th2dt2],'diff(th1,t)=omega1,
+ 'diff(th2,t)=omega2),
+ sol:rkf45(equs,[th1,th2,omega1,omega2],
+ [float(%pi/8),float(%pi/4),0,0],[t,0,8],report=true),
+ plot2d([[discrete,map(lambda([u],part(u,[1,2])),sol)],
+ [discrete,map(lambda([u],part(u,[1,3])),sol)]],
+ [style,[lines,2]],[xlabel,"t (s)"],[ylabel,"angle (rad)"],
+ [legend,"theta1(t)","theta2(t)"])
+)$
+/* ---------------------------------------------------------------------------*/
+( disp("","A mildly stiff problem: The Brusselator."),
+ A:2, B:8.533,
+ sol:rkf45([A+y1^2*y2-(B+1)*y1,B*y1-y1^2*y2],[y1,y2],[1,4.2665],[x,0,20],
+ report=true),
+ plot2d([[discrete,map(lambda([u],part(u,[1,2])),sol)],
+ [discrete,map(lambda([u],part(u,[1,3])),sol)]],[style,[lines,2]],
+ [ylabel,"y1, y2"],[legend,"y1(x)","y2(x)"])
+)$
+/* ---------------------------------------------------------------------------*/
View
185 share/contrib/rkf45/rkf45.mac
@@ -0,0 +1,185 @@
+/*
+--------------------------------------------------------------------------------
+Copyright (C) 2011 Panagiotis Papasotiriou
+
+This software is released under the terms of the GNU General Public License.
+See <http://www.gnu.org/licenses/> for details.
+--------------------------------------------------------------------------------
+ Version: 1
+Date created: September 12, 2011
+ Last update: October 26, 2011
+ Author: Panagiotis J. Papasotiriou
+--------------------------------------------------------------------------------
+Brief description:
+
+rkf45 is a Maxima function for solving initial value problems with automatic
+step size and error control.
+This is an implementation of the Runge-Kutta-Fehlberg 4th-5th-order scheme.
+--------------------------------------------------------------------------------
+Syntax:
+
+rkf45(ode,func,init,interval,options)
+rkf45([ode1,ode2,...],[func1,func2,...],[init1,init2,...],interval,options)
+
+The first form solves a first-order differential equation, ode, with respect to
+the initial condition init, where func is the dependent variable and init is the
+value of the dependent variable at the initial point.
+
+Similarly, the second form solves a system of first-order differential
+equations, ode1,ode2,..., with respect to the initial conditions
+init1,init2,..., where func1,func2,... are the dependent variables and
+init1,init2,... are the values of the dependent variables at the initial point.
+
+Differential equation(s) should be given as expressions depending only on the
+independent and dependent variables, and should define the derivative of the
+dependent variable with respect to the independent variable. For instance, the
+differential equation y'(x)+(x+1)*y=0 should be given as -(x+1)*y.
+The argument "interval" should be a list of three elements. The first element
+identifies the independent variable, while the second and third elements are the
+initial and final values for the independent variable, for instance: [x,0,6].
+Initial value does not need to be less than final value, so an interval like
+[x,6,0] is also valid.
+
+rkf45 accepts the following optional arguments:
+
+* full_solution: A Boolean. If set to true, a full list of the solution at
+ all intermediate points will be returned. If set to false,
+ only the solution at the last integration point is
+ returned. Default: true.
+
+* absolute_tolerance: The desired absolute tolerance of the result. Default:
+ 1e-6.
+
+* max_iterations: Maximum number of iterations. Default: 10000.
+
+* h_start: Initial integration step. Default: one 100th of the
+ integration interval, (interval[3]-interval[2])/100.
+
+* report: A Boolean. If set to true, rkf45 prints a report at
+ exit, giving details about the calculations done.
+ Default: false.
+
+Integration step size ia selected automatically in such a way that the estimated
+local error is less than user-specified absolute tolerance.
+The result is returned as a list with n+1 columns, where n is the number of
+first-order differential equations. The first column contains the values of the
+independent variable selected by the algorithm. The second column contains the
+values of the first dependent variable at the corresponding value of the
+independent variable. Similarly, the third column contains the values of the
+second dependent variable at the corresponding value of the independent
+variable, and so on.
+
+rkf45 can be used to solve moderately stiff initial value problems, although it
+is not designed for that purpose.
+--------------------------------------------------------------------------------
+Examples:
+
+(1) A first-order differential equation, y'=x*(y-1)+3, with y(0)=-2:
+ rkf45(x*(y-1)+3,y,-2,[x,0,4]) returns solution at selected points from x=0
+ to x=4.
+
+(2) A second-order differential equation, y''=x+y*y', with y(-1)=2, y'(-1)=0:
+ rkf45([y2,x+y1*y2],[y1,y2],[2,0],[x,-1,4]) returns solution at selected
+ points from x=-1 to x=4.
+
+See rkf45.dem for more examples. Detailed documentation and several examples
+discussed in detail can be found at the accompanying document rkf45.pdf.
+--------------------------------------------------------------------------------
+*/
+rkf45(odes,funcs,initial,interval,[options]):=block(
+ [numer:true,atol,save_steps,maxit,show_report,xi,xc,yi,h,not_ok,
+ k1,k2,k3,k4,k5,k6,trunc_error,h_optimal,estimated_errors,step_extrema,
+ x_tiny:1e-7*interval[3],bad_steps:0,sol],
+ /* Check interval for errors */
+ if (not(listp(interval))) then error("Error: interval should be a list, but found",interval,"instead."),
+ /* Set optional arguments */
+ atol:assoc('absolute_tolerance,options,1e-6),
+ save_steps:assoc('full_solution,options,true),
+ maxit:assoc('max_iterations,options,10000),
+ show_report:assoc('report,options,false),
+ h:assoc('h_start,options,(interval[3]-interval[2])/100),
+ /* Convert arguments to lists, if necessary */
+ if (not(listp(odes))) then odes:[odes],
+ if (not(listp(funcs))) then funcs:[funcs],
+ if (not(listp(initial))) then initial:[initial],
+ /* Define right-hand-side function */
+ local(f_rhs),
+ define(funmake(f_rhs,cons(interval[1],funcs)),odes),
+ translate(f_rhs),
+ /* Initialize */
+ xi:interval[2], yi:initial,
+ /* Check initial values */
+ if member(false,map(numberp,apply(f_rhs,cons(xi,yi)))) then error("Error: initial should be a list of numbers, but found",funcs,"instead."),
+ /* Set up the rest */
+ estimated_errors:[1e30,-1e30],
+ step_extrema:[1e30,-1e30],
+ if save_steps then sol:[cons(xi,yi)],
+ not_ok:true,
+ /* Main loop */
+ for i:1 thru maxit do (
+ /* Compute solution at xi+h */
+ xc:xi,
+ k1:h*apply(f_rhs,cons(xi,yi)),
+ k2:h*apply(f_rhs,cons(xi+0.25*h,yi+0.25*k1)),
+ k3:h*apply(f_rhs,cons(xi+0.375*h,yi+0.09375*k1+0.28125*k2)),
+ k4:h*apply(f_rhs,cons(xi+9.230769230769231e-1*h,yi+8.793809740555303e-1*k1
+ -3.277196176604461*k2
+ +3.320892125625854*k3)),
+ k5:h*apply(f_rhs,cons(xi+h,yi+2.032407407407407*k1-8*k2+7.173489278752437*k3
+ -2.058966861598441e-1*k4)),
+ k6:h*apply(f_rhs,cons(xi+0.5*h,yi-2.962962962962963e-1*k1+2*k2
+ -1.381676413255361*k3
+ +4.529727095516569e-1*k4-0.275*k5)),
+ /* Estimate local truncation error */
+ trunc_error:lmax(abs(2.777777777777778e-3*k1-2.994152046783626e-2*k3
+ -2.919989367357788e-2*k4+0.02*k5+3.636363636363636e-2*k6
+ ))/abs(h),
+ if (trunc_error<atol) then (
+ /* Step accepted, do it */
+ xi:xc+h,
+ yi:yi+1.157407407407407e-1*k1+5.489278752436647e-1*k3
+ +5.353313840155946e-1*k4-0.2*k5,
+ if save_steps then sol:endcons(cons(xi,yi),sol),
+ estimated_errors[1]:min(trunc_error,estimated_errors[1]),
+ estimated_errors[2]:max(trunc_error,estimated_errors[2]),
+ step_extrema[1]:min(h,step_extrema[1]),
+ step_extrema[2]:max(h,step_extrema[2]),
+ /* Return if final point is reached. */
+ if abs(xi-interval[3])<=x_tiny then (
+ not_ok:false,
+ if not(save_steps) then sol:cons(xi,yi),
+ if show_report then (
+ print("------------------------------------------------------"),
+ print("Info: rkf45:"),
+ print(" Integration points selected:",i+1-bad_steps),
+ print(" Total number of iterations:",i),
+ print(" Bad steps corrected:",bad_steps),
+ print(" Minimum estimated error:",estimated_errors[1]),
+ print(" Maximum estimated error:",estimated_errors[2]),
+ print("Minimum integration step taken:",step_extrema[1]),
+ print("Maximum integration step taken:",step_extrema[2]),
+ print("------------------------------------------------------")
+ ),
+ return()
+ )
+ ) else (
+ /* Step is not accepted, must try again with optimal step */
+ bad_steps:bad_steps+1
+ ),
+ /* Prepare next step */
+ h_optimal:min(max(0.84089641525371*(atol/trunc_error)^0.25,0.1),4)*h,
+ if h_optimal>0 then
+ if xi+h_optimal<interval[3] then h:h_optimal else h:interval[3]-xi
+ else if xi+h_optimal>interval[3] then h:h_optimal else h:interval[3]-xi
+ ),
+ /* Warn user if necessary */
+ if not_ok then (
+ print("Warning: rkf45: Integration stopped at x =",xi,"(stiff problem?)"),
+ print(" Iterations limit has been reached. Check if differential"),
+ print(" equations/initial conditions are given correctly, reduce"),
+ print(" accuracy, and/or increase maximum number of steps.")
+ ),
+ /* Return solution */
+ return(sol)
+)$
+/*----------------------------------------------------------------------------*/
View
BIN share/contrib/rkf45/rkf45.pdf
Binary file not shown.
View
50 share/contrib/rkf45/rtest_rkf45.mac
@@ -0,0 +1,50 @@
+/* ---------------------------------------------------------------------------*/
+/* Test file for rkf45. */
+/* Note: Expected results were obtained on a 64-bit GNU/Linux system. Results
+ obtained on 32-bit systems and/or different platforms may differ slightly.
+ However, float_approx_equal_tolerance is set to 2e-12, which should be
+ enough so that all tests should be passed on any system. */
+/* ---------------------------------------------------------------------------*/
+(kill(all), load('rkf45), float_approx_equal_tolerance:2e-12, 0);
+0$
+/* ---------------------------------------------------------------------------*/
+/* One differential equation of first order. */
+rkf45(-3*x*y^2+1/(x^3+1),y,0,[x,0,5],full_solution=false);
+[5.0,0.039682302083817]$
+/* ---------------------------------------------------------------------------*/
+/* An initial value problem with threshold effect. */
+rkf45(0.206-1.51*f+3.03*f^2/(1+f^2),f,0,[t,0,100],full_solution=false);
+[100.0,1.557094207358533]$
+/* ---------------------------------------------------------------------------*/
+/* A system of two first-order differential equations. */
+rkf45([1.89247311827957*(C-L),52.757793764988*(32-C)+42.20623501199041*(L-C)],
+ [L,C],[150,150],[t,0,2],full_solution=false);
+[2.0,46.84310577179979,38.67012732884611]$
+/* ---------------------------------------------------------------------------*/
+/* A second-order differential equation: the van der Pol equation. */
+rkf45([x2,4*(1-x1^2)*x2-x1],[x1,x2],[0.75,0],[t,0,20],full_solution=false);
+[20.0,1.185928216142604,-0.46868416821173]$
+/* ---------------------------------------------------------------------------*/
+/* A system of two second-order differential equations: the double pendulum. */
+( m1:1, m2:1.5, l1:0.4, l2:0.6, g:9.81,
+ d2th1dt2:(-g*(2*m1+m2)*sin(th1)-m2*g*sin(th1-2*th2)
+ -2*sin(th1-th2)*m2*('diff(th2,t)^2*l2
+ +'diff(th1,t)^2*l1*cos(th1-th2)))
+ /(l1*(2*m1+m2-m2*cos(2*th1-2*th2))),
+ d2th2dt2:(2*sin(th1-th2)*('diff(th1,t)^2*l1*(m1+m2)+g*(m1+m2)*cos(th1)
+ +'diff(th2,t)^2*l2*m2*cos(th1-th2)))
+ /(l2*(2*m1+m2-m2*cos(2*th1-2*th2))),
+ equs:ev([omega1,omega2,d2th1dt2,d2th2dt2],'diff(th1,t)=omega1,
+ 'diff(th2,t)=omega2),
+ sol:rkf45(equs,[th1,th2,omega1,omega2],
+ [float(%pi/8),float(%pi/4),0,0],[t,0,8],full_solution=false)
+);
+[8.0,0.43986621496349,0.31219001091419,-2.057936678736255,-1.092518484540562]$
+/* ---------------------------------------------------------------------------*/
+/* A mildly stiff problem: The Brusselator. */
+rkf45([2+y1^2*y2-(8.533+1)*y1,8.533*y1-y1^2*y2],[y1,y2],[1,4.2665],[x,0,20],
+ full_solution=false);
+[20.0,6.871759807473679,6.850927787700432]$
+/* ---------------------------------------------------------------------------*/
+
+/* batch("/pap/Maxima/rkf45/rtest_rkf45.mac",'test); */
View
4 src/sharefiles.mk
@@ -427,6 +427,10 @@ contrib/rand/vandpol.mac \
contrib/rand/vandpol.usg \
contrib/rand/Wkb \
contrib/README \
+contrib/rkf45/rkf45.dem \
+contrib/rkf45/rkf45.mac \
+contrib/rkf45/rkf45.pdf \
+contrib/rkf45/rtest_rkf45.mac \
contrib/sarag/aliases.mac \
contrib/sarag/arag_test.mac \
contrib/sarag/certificateOfPositivity.mac \

0 comments on commit c2c9482

Please sign in to comment.
Something went wrong with that request. Please try again.