Differential equations (Runge-Kutta-Gill)

(Back)
Email the author:Paolo Bonzini

Usage:

- Edit the program after label EQ. Enter a statement like this: {<equation for y'>,<equation for y''>,...)=>L2 The n-th derivative of y is written as L1(n+1) and y is written as L1(1). If you have no equation for the n-th derivative, write L1(n+1) as the n-th element of L2.

For example y''=y'*y'/y+y becomes {L1(2),L1(2)^2/L1(1)+L1(1)}=>L2

- Then run the program and answer its questions. mat A, at the end, contains the results in this form: ((X,Y,Y',...), (X,Y,Y',...),...)

The program (\ is square root):

Goto START
Label EQ
{L1(2),L1(2)^2/L1(1)+L1(1)}=>L2
Return

Label START
Print "STEP"
Input B
Print "START X"
Input X
Print "START Y"
Input L1
Print "END X"
Input E
dim(L1)=>dim(L3)
list->mat(augment({X},L1),mat A)
Gosub EQ

Label LOOP
BL2=>L4
L1+.5L4-L3=>L1
L4-2L3=>L3

X+.5B=>X
1-\.5=>V
Gosub EQ
BL2=>L4
(L4-L3)V=>L5
L1+L5=>L1
L3+3L5-VL4=>L3

1+\.5=>V
BL2=>L4
(L4-L3)V=>L5
L1+L5=>L1
L3+3L5-VL4=>L3

X+.5B=>X
Gosub EQ
BL2=>L4
(L4-2L3)/6=>L5
L1+L5=>L1
L3+3L5-.5L4=>L3

Print X
list->mat(augment({X},L1),mat B)
augment(mat A,mat B)=>mat A
If X<E Goto LOOP
trans mat A=>mat A