/* -*-MACSYMA-*- Simple programs for the numerical solution of ordinary differential equations. A better implementation of this stuff is NDIFFQ. -GJC */ /* The main use for these RUNGE-KUTTA formulas is to get starting values for predictor-corrector methods. */ /* 11:57am Tuesday, 8 July 1980 George Carrette.*/ /* This is meant to be a useful program, but even more so it is meant to show how to use modedeclarations and compilation to write numerical programs. TRANSLATE_FILE("SHARE2\;DIFFEQ"), and then COMPILE_LISP_FILE("SHARE2\;DIFFEQ"). */ EVAL_WHEN(translate,transcompile:true /* ,packagefile:true */)$ /* HERALD_PACKAGE(DIFFEQ_RUNGE)$ */ FLOATCHECK(X):=(MODEDECLARE(X,ANY,function(float),float), X:FLOAT(X), IF NOT NUMBERP(X) THEN ERROR(X,"not a floating point number."), X)$ MODEDECLARE(FUNCTION(FLOATCHECK),FLOAT)$ DEFINE_VARIABLE(X_WE_ARE_CALCULATING,"WE ARE NOT CALCULATING",ANY)$ DEFINE_VARIABLE(DIFFEQ_RUNTIME,apply('STATUS,'[RUNTIME]),ANY)$ /* can't use control-U on VAX DIFFEQ_TTYINT_FUN():= IF X_WE_ARE_CALCULATING="WE ARE NOT CALCULATING" THEN FALSE ELSE PRINT(" Calculating at X=",X_WE_ARE_CALCULATING, apply('STATUS,'[RUNTIME])-DIFFEQ_RUNTIME,"so far.")$ */ EVAL_WHEN([TRANSLATE,BATCH,DEMO], /* a macro to set up a bunch of calls to FLOATCHECK. */ FLOATCHECKL([L])::= FUNMAKE("(", MAPLIST(LAMBDA([VAR], FUNMAKE(":",[VAR,FUNMAKE('FLOATCHECK,[VAR])])), L)) /* , */ /* A macro to bind the ttyintfun */ /* WITH_DIFFEQ_TTYINTFUN([L])::= BUILDQ([L], BLOCK([ttyintfun:'diffeq_ttyint_fun,diffeq_runtime: apply('status,'[runtime]), X_WE_ARE_CALCULATING:"WE ARE NOT CALCULATING"], SPLICE(L))) */ )$ /* RUNGE1(F,X0,X1,DX,Y0) returns a list of two lists the X points and the Y points. The following formuli are used: dY -- = F(X, Y) dX K K K K 4 3 2 1 Y = Y + -- + -- + -- + -- N + 1 N 6 3 3 6 K = H F(X , Y ) 1 N N K H 1 K = H F(X + -, Y + --) 2 N 2 N 2 K H 2 K = H F(X + -, Y + --) 3 N 2 N 2 K = H F(X + H, Y + K ) 4 N N 3 */ RUNGE1(F,X0,X1,H,Y0):= (MODEDECLARE([F,X0,X1,H,Y0],ANY), /* This user interface function provides error checking. When compiled, the function RUNGE1_INTERNAL cannot check that the type of its arguments is float because of the flonum declaration it uses. */ FLOATCHECKL(X0,X1,H,Y0), /* WITH_DIFFEQ_TTYINTFUN( */ RUNGE1_INTERNAL(F,X0,X1,H,Y0) /* ) */ )$ RUNGE1_INTERNAL(F,X0,X1,H,Y0):= (MODEDECLARE([FUNCTION(F),X0,X1,H,Y0],FLOAT,F,ANY), BLOCK([Y_LIST:[],X_LIST:[],YP_LIST:[],K1,K2,K3,K4], MODEDECLARE([Y_LIST,X_LIST,YP_LIST],LIST,[X,K1,K2,K3,K4],FLOAT), FOR X:X0 THRU X1 STEP H DO (X_WE_ARE_CALCULATING:X, Y_LIST:CONS(Y0,Y_LIST), X_LIST:CONS(X,X_LIST), K1:APPLY(F,[X,Y0]), /* using K1 as a temp register. */ YP_LIST:CONS(K1,YP_LIST), K1:H*K1, /* the real K1. */ K2:H*APPLY(F,[X+H/2,Y0+K1/2]), K3:H*APPLY(F,[X+H/2,Y0+K2/2]), K4:H*APPLY(F,[X+H,Y0+K3]), Y0:Y0+(K1+K2)/6+(K2+K3)/3), X_WE_ARE_CALCULATING:"WE ARE NOT CALCULATING", [REVERSE(X_LIST),REVERSE(Y_LIST),REVERSE(YP_LIST)]))$ /* RUNGE2(F,X0,X1,DX,Y0,Y_PRIME0), returns a list [ , , ] ref: abromowitz & stegun pp 897. 2 d Y dY --- = F(X, Y, --) 2 dX dX K + K + K dY 3 2 1 Y = H ((--) + ------------) + Y N + 1 dX N 6 N K + 2 K + 2 K + K dY dY 4 3 2 1 (--) = (--) + ---------------------- dX N + 1 dX N 6 dY K = H F(X , Y , (--) ) 1 N N dX N dY H (--) K H K H dX N 1 dY 1 K = H F(X + -, ------- + Y + ----, (--) + --) 2 N 2 2 N 8 dX N 2 dY H (--) K H K H dX N 1 dY 2 K = H F(X + -, ------- + Y + ----, (--) + --) 3 N 2 2 N 8 dX N 2 K H dY 3 dY K = H F(X + H, H (--) + Y + ----, (--) + K ) 4 N dX N N 2 dX N 3 */ RUNGE2(F,X0,X1,H,Y0,YP0):= (MODEDECLARE([F,X0,X1,H,Y0,YP0],ANY), FLOATCHECKL(X0,X1,H,Y0,YP0), /* WITH_DIFFEQ_TTYINTFUN( */ RUNGE2_INTERNAL(F,X0,X1,H,Y0,YP0) /* ) */)$ RUNGE2_INTERNAL(F,X0,X1,H,Y0,YP0):= (MODEDECLARE([FUNCTION(F),X0,X1,H,Y0,YP0],FLOAT,F,ANY), BLOCK([Y_LIST:[],X_LIST:[],YP_LIST:[],YPP_LIST:[],K1,K2,K3,K4,TEMP], MODEDECLARE([Y_LIST,X_LIST,YP_LIST,YPP_LIST],LIST, [X,K1,K2,K3,K4,TEMP],FLOAT), FOR X:X0 THRU X1 STEP H DO (X_WE_ARE_CALCULATING:X, Y_LIST:CONS(Y0,Y_LIST), X_LIST:CONS(X,X_LIST), YP_LIST:CONS(YP0,YP_LIST), K1:APPLY(F,[X ,Y0 ,YP0 ]), YPP_LIST:CONS(K1,YPP_LIST), K1:H*K1, K2:H*APPLY(F,[X+H/2,Y0+H/2*(YP0+K1/4),YP0+K1/2]), K3:H*APPLY(F,[X+H/2,Y0+H/2*(YP0+K2/4),YP0+K2/2]), K4:H*APPLY(F,[X+H ,Y0+H *(YP0+K3/2),YP0+K3 ]), /* This next must be a parallel assignment. */ TEMP:YP0+(K1+2*(K2+K3)+K4)/6, Y0:Y0+H*(YP0+(K1+K2+K3)/6), YP0:TEMP), X_WE_ARE_CALCULATING:"WE ARE NOT CALCULATING", [REVERSE(X_LIST),REVERSE(Y_LIST),REVERSE(YP_LIST),REVERSE(YPP_LIST)]))$ /* This handles N first order equations RUNGEN([F1,F2,F3],XA,XB,H,[Y1A,Y2A,Y3A]) The functions, e.g. DY1/DX=F1(X,Y1,Y2,Y3) are functions of N+1 arguments. */ RUNGEN(FL,XA,XB,H,YAL):= (MODEDECLARE([FL,XA,XB,H,YAL],any), BLOCK([ORDER:LENGTH(FL)], MODEDECLARE(ORDER,FIXNUM), IF ORDER#LENGTH(YAL) THEN ERROR("Wrong number of initial values",FL,YAL), FLOATCHECKL(XA,XB,H), YAL:MAPLIST('FLOATCHECK,YAL), /* WITH_DIFFEQ_TTYINTFUN( */ RUNGEN_INTERNAL(FL,XA,XB,H,YAL) /* ) */ ))$ /* VAPPLY is an inefficient but general way to solve this programming problem. It is inefficient because macsyma lists make inefficient vectors.*/ VAPPLY(FL,ARGL):= (MODEDECLARE([FL,ARGL],LIST), DECLARE(ARGL,SPECIAL), MAPLIST(LAMBDA([F],MODEDECLARE(F,ANY),APPLY(F,ARGL)),FL))$ DECLARE(LISTARITH,SPECIAL)$ RUNGEN_INTERNAL(F,XA,XB,H,YA):= (MODEDECLARE([XA,XB,H],FLOAT,[F,YA],LIST), BLOCK([Y_LIST:[],X_LIST:[],YP_LIST:[],K1,K2,K3,K4,LISTARITH:TRUE], MODEDECLARE([Y_LIST,X_LIST,YP_LIST,K1,K2,K3,K4],LIST, [X],FLOAT), FOR X:XA THRU XB STEP H DO (X_WE_ARE_CALCULATING:X, Y_LIST:CONS(YA,Y_LIST), X_LIST:CONS(X,X_LIST), K1:VAPPLY(F,CONS(X,YA)), /* using K1 as a temp register. */ YP_LIST:CONS(K1,YP_LIST), K1:H*K1, /* the real K1. */ K2:H*VAPPLY(F,CONS(X+H/2,YA+K1/2.0)), K3:H*VAPPLY(F,CONS(X+H/2,YA+K2/2.0)), K4:H*VAPPLY(F,CONS(X+H,YA+K3)), YA:YA+(K1+K2)/6.0+(K2+K3)/3.0), X_WE_ARE_CALCULATING:"WE ARE NOT CALCULATING", [REVERSE(X_LIST),REVERSE(Y_LIST),REVERSE(YP_LIST)]))$