share/numeric/diffeq.mac

   1 /* -*-MACSYMA-*-
   2    Simple programs for the numerical solution of ordinary
   3    differential equations.
   4
   5    A better implementation of this stuff is NDIFFQ. -GJC
   6
   7  */
   8
   9
  10 /* The main use for these RUNGE-KUTTA formulas is to get starting
  11    values for predictor-corrector methods. */
  12
  13 /* 11:57am  Tuesday, 8 July 1980 George Carrette.*/
  14
  15 /* This is meant to be a useful program, but even more so it
  16    is meant to show how to use modedeclarations and compilation
  17    to write numerical programs.
  18
  19    TRANSLATE_FILE("SHARE2\;DIFFEQ"), and then
  20    COMPILE_LISP_FILE("SHARE2\;DIFFEQ").
  21 */
  22
  23 status([x]):='done$
  24
  25 /*eval_when(translate,packagefile:true)$*/
  26
  27 /* HERALD_PACKAGE(DIFFEQ_RUNGE)$ */
  28
  29 floatcheck(x):=(modedeclare(x,any,function(float),float),
  30    x:float(x),
  31    if not numberp(x) then error(x,"not a floating point number."),
  32    x)$
  33
  34 modedeclare(function(floatcheck),float)$
  35
  36 define_variable(x_we_are_calculating,"we are not calculating",any)$
  37 define_variable(diffeq_runtime,apply('status,'[runtime]),any)$
  38
  39 /* can't use control-U on VAX
  40 DIFFEQ_TTYINT_FUN():=
  41  IF X_WE_ARE_CALCULATING="WE ARE NOT CALCULATING"
  42    THEN FALSE
  43    ELSE PRINT("
  44 Calculating at X=",X_WE_ARE_CALCULATING,
  45         apply('STATUS,'[RUNTIME])-DIFFEQ_RUNTIME,"so far.")$ */
  46
  47
  48 eval_when([translate,batch,demo],
  49    /* a macro to set up a bunch of calls to FLOATCHECK. */
  50    floatcheckl([l])::=
  51     funmake("(",
  52      maplist(lambda([var],
  53                funmake(":",[var,funmake('floatcheck,[var])])),
  54              l)) /* , */
  55    /* A macro to bind the ttyintfun */
  56 /*   WITH_DIFFEQ_TTYINTFUN([L])::=
  57     BUILDQ([L],
  58       BLOCK([ttyintfun:'diffeq_ttyint_fun,diffeq_runtime:
  59 apply('status,'[runtime]),
  60              X_WE_ARE_CALCULATING:"WE ARE NOT CALCULATING"],
  61             SPLICE(L))) */
  62 )$
  63
  64 /* RUNGE1(F,X0,X1,DX,Y0) returns a list of two lists
  65    the X points and the Y points. The following formuli
  66    are used:
  67
  68            dY
  69            -- = F(X, Y)
  70            dX
  71
  72                 K    K    K    K
  73                  4    3    2    1
  74   Y      = Y  + -- + -- + -- + --
  75    N + 1    N   6    3    3    6
  76
  77
  78    K  = H F(X , Y )
  79     1        N   N
  80
  81                          K
  82                  H        1
  83    K  = H F(X  + -, Y  + --)
  84     2        N   2   N   2
  85
  86                          K
  87                  H        2
  88    K  = H F(X  + -, Y  + --)
  89     3        N   2   N   2
  90
  91
  92    K  = H F(X  + H, Y  + K )
  93     4        N       N    3
  94 */
  95
  96 runge1(f,x0,x1,h,y0):=
  97 (modedeclare([f,x0,x1,h,y0],any),
  98  /* This user interface function provides error checking.
  99     When compiled, the function RUNGE1_INTERNAL cannot check that the
 100     type of its arguments is float because of the flonum declaration
 101     it uses. */
 102   floatcheckl(x0,x1,h,y0),
 103 /*  with_diffeq_ttyintfun( */ runge1_internal(f,x0,x1,h,y0) /* ) */ )$
 104
 105 runge1_internal(f,x0,x1,h,y0):=
 106  (modedeclare([function(f),x0,x1,h,y0],float,f,any),
 107   block([y_list:[],x_list:[],yp_list:[],k1,k2,k3,k4],
 108         modedeclare([y_list,x_list,yp_list],list,[x,k1,k2,k3,k4],float),
 109         for x:x0 thru x1 step h
 110          do (x_we_are_calculating:x,
 111              y_list:cons(y0,y_list),
 112              x_list:cons(x,x_list),
 113              k1:apply(f,[x,y0]),         /* using k1 as a temp register. */
 114              yp_list:cons(k1,yp_list),
 115              k1:h*k1,                    /* the real k1. */
 116              k2:h*apply(f,[x+h/2,y0+k1/2]),
 117              k3:h*apply(f,[x+h/2,y0+k2/2]),
 118              k4:h*apply(f,[x+h,y0+k3]),
 119              y0:y0+(k1+k2)/6+(k2+k3)/3),
 120         x_we_are_calculating:"we are not calculating",
 121         [reverse(x_list),reverse(y_list),reverse(yp_list)]))$
 122
 123 /* RUNGE2(F,X0,X1,DX,Y0,Y_PRIME0), returns a list [ <x's> , <y's> , <dy/dx's> ]
 124 ref: abromowitz & stegun pp 897.
 125
 126     2
 127  d Y           dY
 128  --- = F(X, Y, --)
 129    2           dX
 130  dX
 131
 132
 133                        K  + K  + K
 134                dY       3    2    1
 135   Y      = H ((--)  + ------------) + Y
 136    N + 1       dX N        6           N
 137
 138                        K  + 2 K  + 2 K  + K
 139    dY          dY       4      3      2    1
 140   (--)      = (--)  + ----------------------
 141    dX N + 1    dX N             6
 142
 143
 144                     dY
 145   K  = H F(X , Y , (--) )
 146    1        N   N   dX N
 147
 148
 149                      dY
 150                   H (--)         K  H          K
 151                H     dX N         1     dY      1
 152  K  = H F(X  + -, ------- + Y  + ----, (--)  + --)
 153   2        N   2     2       N    8     dX N   2
 154
 155
 156                       dY
 157                    H (--)        K  H           K
 158                 H     dX N        1      dY      2
 159   K  = H F(X  + -, ------- + Y  + ----, (--)  + --)
 160    3        N   2     2      N    8     dX N   2
 161
 162                                   K  H
 163                       dY           3     dY
 164   K  = H F(X  + H, H (--)  + Y  + ----, (--)  + K )
 165    4        N         dX N    N    2     dX N    3
 166
 167
 168 */
 169
 170 runge2(f,x0,x1,h,y0,yp0):=
 171  (modedeclare([f,x0,x1,h,y0,yp0],any),
 172   floatcheckl(x0,x1,h,y0,yp0),
 173 /*  with_diffeq_ttyintfun( */ runge2_internal(f,x0,x1,h,y0,yp0) /* ) */)$
 174
 175 runge2_internal(f,x0,x1,h,y0,yp0):=
 176  (modedeclare([function(f),x0,x1,h,y0,yp0],float,f,any),
 177   block([y_list:[],x_list:[],yp_list:[],ypp_list:[],k1,k2,k3,k4,temp],
 178         modedeclare([y_list,x_list,yp_list,ypp_list],list,
 179                     [x,k1,k2,k3,k4,temp],float),
 180         for x:x0 thru x1 step h
 181          do (x_we_are_calculating:x,
 182              y_list:cons(y0,y_list),
 183              x_list:cons(x,x_list),
 184              yp_list:cons(yp0,yp_list),
 185              k1:apply(f,[x    ,y0               ,yp0     ]),
 186              ypp_list:cons(k1,ypp_list),
 187              k1:h*k1,
 188              k2:h*apply(f,[x+h/2,y0+h/2*(yp0+k1/4),yp0+k1/2]),
 189              k3:h*apply(f,[x+h/2,y0+h/2*(yp0+k2/4),yp0+k2/2]),
 190              k4:h*apply(f,[x+h  ,y0+h  *(yp0+k3/2),yp0+k3  ]),
 191              /* this next must be a parallel assignment. */
 192              temp:yp0+(k1+2*(k2+k3)+k4)/6,
 193              y0:y0+h*(yp0+(k1+k2+k3)/6),
 194              yp0:temp),
 195         x_we_are_calculating:"we are not calculating",
 196         [reverse(x_list),reverse(y_list),reverse(yp_list),reverse(ypp_list)]))$
 197 \f
 198 /* This handles N first order equations
 199    RUNGEN([F1,F2,F3],XA,XB,H,[Y1A,Y2A,Y3A])
 200    The functions, e.g. DY1/DX=F1(X,Y1,Y2,Y3) are functions of N+1 arguments.
 201 */
 202
 203 rungen(fl,xa,xb,h,yal):=
 204  (modedeclare([fl,xa,xb,h,yal],any),
 205   block([order:length(fl)],
 206         modedeclare(order,fixnum),
 207         if order#length(yal) then error("wrong number of initial values",fl,yal),
 208         floatcheckl(xa,xb,h),
 209         yal:maplist('floatcheck,yal),
 210 /*      with_diffeq_ttyintfun( */ rungen_internal(fl,xa,xb,h,yal) /* ) */ ))$
 211
 212 /* VAPPLY is an inefficient but general way to solve this programming
 213    problem. It is inefficient because macsyma lists make inefficient
 214    vectors.*/
 215
 216 vapply(fl,argl):=
 217  (modedeclare([fl,argl],list),
 218   maplist(lambda([f],modedeclare(f,any),apply(f,argl)),fl))$
 219
 220 rungen_internal(f,xa,xb,h,ya):=
 221  (modedeclare([xa,xb,h],float,[f,ya],list),
 222   block([y_list:[],x_list:[],yp_list:[],k1,k2,k3,k4,listarith:true],
 223         modedeclare([y_list,x_list,yp_list,k1,k2,k3,k4],list,
 224                     [x],float),
 225         for x:xa thru xb step h
 226          do (x_we_are_calculating:x,
 227              y_list:cons(ya,y_list),
 228              x_list:cons(x,x_list),
 229              k1:vapply(f,cons(x,ya)),         /* using k1 as a temp register. */
 230              yp_list:cons(k1,yp_list),
 231              k1:h*k1,                    /* the real k1. */
 232              k2:h*vapply(f,cons(x+h/2,ya+k1/2.0)),
 233              k3:h*vapply(f,cons(x+h/2,ya+k2/2.0)),
 234              k4:h*vapply(f,cons(x+h,ya+k3)),
 235              ya:ya+(k1+k2)/6.0+(k2+k3)/3.0),
 236         x_we_are_calculating:"we are not calculating",
 237         [reverse(x_list),reverse(y_list),reverse(yp_list)]))$
 238
 239
 240