share/contrib/rkf45/rtest_rkf45.mac

   1 /* ---------------------------------------------------------------------------*/
   2 /* Test file for rkf45.                                                       */
   3 /* Note: Expected results were obtained on a 64-bit GNU/Linux system. Results
   4    obtained on 32-bit systems and/or different platforms may differ slightly.
   5    However, float_approx_equal_tolerance is set to 2e-12, which should be
   6    enough so that all tests should be passed on any system.                   */
   7 /* ---------------------------------------------------------------------------*/
   8 (kill(all), load('rkf45), float_approx_equal_tolerance:2e-12, 0);
   9 0$
  10 /* ---------------------------------------------------------------------------*/
  11 /* One differential equation of first order.                                  */
  12 rkf45(-3*x*y^2+1/(x^3+1),y,0,[x,0,5],full_solution=false);
  13 [5.0,0.039682302083817]$
  14 /* ---------------------------------------------------------------------------*/
  15 /* An initial value problem with threshold effect.                            */
  16 rkf45(0.206-1.51*f+3.03*f^2/(1+f^2),f,0,[t,0,100],full_solution=false);
  17 [100.0,1.557094207358533]$
  18 /* ---------------------------------------------------------------------------*/
  19 /* A system of two first-order differential equations.                        */
  20 rkf45([1.89247311827957*(C-L),52.757793764988*(32-C)+42.20623501199041*(L-C)],
  21       [L,C],[150,150],[t,0,2],full_solution=false);
  22 [2.0,46.84310577179979,38.67012732884611]$
  23 /* ---------------------------------------------------------------------------*/
  24 /* A second-order differential equation: the van der Pol equation.            */
  25 rkf45([x2,4*(1-x1^2)*x2-x1],[x1,x2],[0.75,0],[t,0,20],full_solution=false);
  26 [20.0,1.185928216142604,-0.46868416821173]$
  27 /* ---------------------------------------------------------------------------*/
  28 /* A system of two second-order differential equations: the double pendulum.  */
  29 ( m1:1, m2:1.5, l1:0.4, l2:0.6, g:9.81,
  30   d2th1dt2:(-g*(2*m1+m2)*sin(th1)-m2*g*sin(th1-2*th2)
  31             -2*sin(th1-th2)*m2*('diff(th2,t)^2*l2
  32             +'diff(th1,t)^2*l1*cos(th1-th2)))
  33            /(l1*(2*m1+m2-m2*cos(2*th1-2*th2))),
  34   d2th2dt2:(2*sin(th1-th2)*('diff(th1,t)^2*l1*(m1+m2)+g*(m1+m2)*cos(th1)
  35             +'diff(th2,t)^2*l2*m2*cos(th1-th2)))
  36            /(l2*(2*m1+m2-m2*cos(2*th1-2*th2))),
  37   equs:ev([omega1,omega2,d2th1dt2,d2th2dt2],'diff(th1,t)=omega1,
  38           'diff(th2,t)=omega2),
  39   sol:rkf45(equs,[th1,th2,omega1,omega2],
  40             [float(%pi/8),float(%pi/4),0,0],[t,0,8],full_solution=false)
  41 );
  42 [8.0,0.43986621496349,0.31219001091419,-2.057936678736255,-1.092518484540562]$
  43 /* ---------------------------------------------------------------------------*/
  44 /* A mildly stiff problem: The Brusselator.                                   */
  45 rkf45([2+y1^2*y2-(8.533+1)*y1,8.533*y1-y1^2*y2],[y1,y2],[1,4.2665],[x,0,20],
  46       full_solution=false);
  47 [20.0,6.871759807473679,6.850927787700432]$
  48 /* ---------------------------------------------------------------------------*/
  49
  50 /* see thread http://permalink.gmane.org/gmane.comp.mathematics.maxima.general/45146
  51    -constants in ode need to be coerced to floats
  52    -initial and interval need to be coerced to floats
  53 */
  54 block([vx,x,t,dx_dt,dvx_dt,ode],
  55   dx_dt : vx,
  56   dvx_dt : -4*%pi^2*x + 100*sin(4*%pi*t + %pi/4)/(0.5 + t^2),
  57   ode : [dx_dt,  dvx_dt],
  58   errcatch(rkfsoln : rkf45 (ode, [x, vx], [1, 0], [t, 0, 1])),
  59   is(%% = []));
  60 false $
  61
  62 /* Error handling for improper interval argument */
  63 errcatch (rkf45(y,y,1,[t,0]));
  64 []$
  65 errcatch (rkf45(y,y,1,t,0,1));
  66 []$
  67
  68 /* see thread http://permalink.gmane.org/gmane.comp.mathematics.maxima.general/45195
  69 */
  70 is( abs(last(last(rkf45 (-t*y,y,1,[t,0,1,0.01]))) - exp(-0.5)) < 1e-6 );
  71 true $
  72
  73 (reset(float_approx_equal_tolerance),0);
  74 0$
  75
  76 /* batch("/pap/Maxima/rkf45/rtest_rkf45.mac",'test); */