share/numeric/diffeq.dem

   1 /* a demo of solving numerical differential equations and
   2    using PLOT2 to plot the results.
   3
   4    N.B. A better implementation is in NDIFFQ -GJC
   5
   6    This demo must be run in a BARE BARE fresh macsyma.
   7
   8    This also illustrates the use of compilation and MODEDECLARATION.
   9
  10    1:24pm  Tuesday, 8 July 1980 George Joseph Carrette.
  11 */
  12 /*
  13 (DYNAMALLOC:TRUE,LOADFILE(NUMER,AUTOLOAD,DSK,NUMER))$
  14
  15 LOAD_PACKAGE(DIFFEQ,"SHARE2\;DIFFEQ")$
  16 */
  17 load("diffeq")$
  18 if showtime=false then showtime:'all$
  19
  20 /* The van der Pol equation. */
  21
  22 'diff('u,'t,2)='micro*(1-'u^2)*'diff('u,'t)-'u;
  23
  24 /* RUNGE2(F,X0,X1,H,Y0,YP0) is the call for a second order equation. */
  25
  26 define_variable(micro,1.0,float)$
  27
  28 vander(t,y,yp):=(mode_declare([t,y,yp],float),micro*(1-y^2)*yp-y);
  29 /*
  30 COMPILE("DSK:NUMER\;.TEMP. DEMO",VANDER)$
  31 */
  32 /* type VANDER and a semi-colon */
  33 compile();
  34 /* It is interesting to plot a few paths with different values of
  35    micro. */
  36 vanderpol&&
  37
  38 /* You can type control-U during the calculation to see how it
  39    is progressing. */
  40
  41 (x0:0.0,x1:10.0,npoints:100,dx:(x1-x0)/npoints)$
  42
  43 vpm1:(micro:0.1,runge2('vander,x0,x1,dx,0.1,0.1))$
  44 (vpx:vpm1[1],vpm1:rest(vpm1))$
  45 vpm2:rest((micro:0.5,runge2('vander,x0,x1,dx,0.1,0.1)))$
  46 vpm3:rest((micro:1.5,runge2('vander,x0,x1,dx,0.1,0.1)))$
  47 vpm4:rest((micro:3.5,runge2('vander,x0,x1,dx,0.1,0.1)))$
  48
  49 /* the plot will be phase space, Y vs. Y' */
  50 /* not in DOE MACSYMA at the present time
  51  AP(X,Y):=APPLY('GRAPH2,[X,Y,[0,1,3,8]])$ */
  52
  53 ap(x,y):=(apply('graph,[x[1],y[1],["*"]]),
  54 apply('graph,[x[2],y[2],["+"]]),
  55 apply('graph,[x[3],y[3],["^"]]),
  56 apply('graph,[x[4],y[4],["-"]]))$
  57
  58 ap([vpm1[1],vpm2[1],vpm3[1],vpm4[1]],
  59    [vpm1[2],vpm2[2],vpm3[2],vpm4[2]]);
  60
  61
  62         /* A time X vs. Y graph. */
  63
  64 ap([vpx,vpx,vpx,vpx],[vpm1[1],vpm2[1],vpm3[1],vpm4[1]]);
  65
  66
  67         /* coupled equations. */
  68 kill(arrays,values,functions,labels)$
  69 spring&&
  70
  71 'diff('x,'t,2)=-'k/'m*'x;
  72
  73 /* reducing this to a system of first order equations ...
  74    (somebody should write a macsyma program to do this
  75     automatically!) */
  76
  77 define_variable(k,1.0,float)$
  78 define_variable(m,1.0,float)$
  79
  80 dvel(t,x,vel):=(modedeclare([t,x,vel],float),-k/m*x)$
  81 dx(t,x,vel):=(modedeclare([t,x,vel],float),vel)$
  82 /*
  83 COMPILE("DSK:NUMER\;.TEMP. DEMO",DVEL,DX)$
  84 */
  85 /* answer '[DVEL,DX] followed by a semi-colon */
  86 compile();
  87 /* The general RUNGEN is not very efficient, it conses up
  88    a lot of temporary lists to serve as vectors. So this
  89    simple problem will take about 10 seconds CPU time.
  90    Type control-U during the calculation to see how
  91    far it has gotten. */
  92
  93 sol:rungen('[dx,dvel],0,10,.15,'[1,0])$
  94
  95 /* The solution is of the form
  96    [ <list-of-t> , <list of [X,VEL] pairs>, <list of [X',VEL']> ]
  97    of course in this case the X' is the same as VEL. The
  98    generality leads to a little waste in this case. RUNGE2 is
  99    better for solving second order equations. */
 100
 101 t:sol[1]$
 102 x:maplist(lambda([u],u[1]),sol[2])$
 103 vel:maplist(lambda([u],u[2]),sol[2])$
 104
 105 /* with K=1 and M=1 we should see a sine and cosine wave of
 106    period 2*%PI, since I picked the initial conditions X=1 and VEL=0
 107 */
 108
 109 graph(t,[x,vel],["*","+"])$
 110
 111 /* end of demo */
 112 "end of demo"$