share/contrib/rkf45/rkf45.dem

   1 /* ---------------------------------------------------------------------------*/
   2 /* Demo file for rkf45.                                                       */
   3 /* ---------------------------------------------------------------------------*/
   4 load('rkf45)$
   5 /* ---------------------------------------------------------------------------*/
   6 ( disp("","One differential equation of first order."),
   7   sol:rkf45(-3*x*y^2+1/(x^3+1),y,0,[x,0,5],report=true),
   8   plot2d([discrete,sol],[style,[lines,2]])
   9 )$
  10 /* ---------------------------------------------------------------------------*/
  11 ( disp("","An initial value problem with threshold effect."),
  12   sol:makelist(rkf45(equ,f,0,[t,0,100]),equ,
  13                makelist(s-1.51*f+3.03*f^2/(1+f^2),s,[0.206,0.204,0.202])),
  14   plot2d(makelist([discrete,s],s,sol),[style,[lines,2]],[xlabel,"t"],
  15          [ylabel,"f"],[legend,"s=0.206","s=0.204","s=0.202"],
  16          [gnuplot_preamble,"set key left"])
  17 )$
  18 /* ---------------------------------------------------------------------------*/
  19 ( disp("","A system of two first-order differential equations."),
  20   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,
  21   sol:rkf45([k1*(C-L),k2*(32-C)+k3*(L-C)],[L,C],[150,150],[t,0,2],report=true),
  22   plot2d([[discrete,map(lambda([u],part(u,[1,2])),sol)],
  23           [discrete,map(lambda([u],part(u,[1,3])),sol)]],
  24          [style,[lines,2]],[xlabel,"time (hours)"],[ylabel,"temperature (F)"],
  25          [legend,"liquid, L(t)","container, C(t)"])
  26 )$
  27 /* ---------------------------------------------------------------------------*/
  28 ( disp("","A second-order differential equation: the van der Pol equation."),
  29   equ:'diff(x,t,2)+4*(x^2-1)*'diff(x,t)+x=0,
  30   equ2:['diff(x1,t)=x2,ev(solve(equ,'diff(x,t,2))[1],'diff(x,t,2)='diff(x2,t),
  31                                                      'diff(x,t)=x2,x=x1)],
  32   equ2:map(rhs,equ2),
  33   sol:rkf45(equ2,[x1,x2],[0.75,0],[t,0,20],report=true),
  34   plot2d([[discrete,map(lambda([u],part(u,[1,2])),sol)],
  35           [discrete,map(lambda([u],part(u,[1,3])),sol)]],
  36          [style,[lines,2]],[xlabel,"t"],[ylabel,"x, x'"],
  37          [legend,"x(t)","x'(t)"])
  38 )$
  39 /* ---------------------------------------------------------------------------*/
  40 ( disp("","A system of two second-order differential equations:",
  41        "the double pendulum."),
  42   m1:1, m2:1.5, l1:0.4, l2:0.6, g:9.81,
  43   d2th1dt2:(-g*(2*m1+m2)*sin(th1)-m2*g*sin(th1-2*th2)
  44             -2*sin(th1-th2)*m2*('diff(th2,t)^2*l2
  45             +'diff(th1,t)^2*l1*cos(th1-th2)))
  46            /(l1*(2*m1+m2-m2*cos(2*th1-2*th2))),
  47   d2th2dt2:(2*sin(th1-th2)*('diff(th1,t)^2*l1*(m1+m2)+g*(m1+m2)*cos(th1)
  48             +'diff(th2,t)^2*l2*m2*cos(th1-th2)))
  49            /(l2*(2*m1+m2-m2*cos(2*th1-2*th2))),
  50   equs:ev([omega1,omega2,d2th1dt2,d2th2dt2],'diff(th1,t)=omega1,
  51           'diff(th2,t)=omega2),
  52   sol:rkf45(equs,[th1,th2,omega1,omega2],
  53             [float(%pi/8),float(%pi/4),0,0],[t,0,8],report=true),
  54   plot2d([[discrete,map(lambda([u],part(u,[1,2])),sol)],
  55           [discrete,map(lambda([u],part(u,[1,3])),sol)]],
  56          [style,[lines,2]],[xlabel,"t (s)"],[ylabel,"angle (rad)"],
  57          [legend,"theta1(t)","theta2(t)"])
  58 )$
  59 /* ---------------------------------------------------------------------------*/
  60 ( disp("","A mildly stiff problem: The Brusselator."),
  61   A:2, B:8.533,
  62   sol:rkf45([A+y1^2*y2-(B+1)*y1,B*y1-y1^2*y2],[y1,y2],[1,4.2665],[x,0,20],
  63       report=true),
  64   plot2d([[discrete,map(lambda([u],part(u,[1,2])),sol)],
  65           [discrete,map(lambda([u],part(u,[1,3])),sol)]],[style,[lines,2]],
  66           [ylabel,"y1, y2"],[legend,"y1(x)","y2(x)"])
  67 )$
  68 /* ---------------------------------------------------------------------------*/