Numerically solving Ordinary Differential Equations
My first attempt at an interactive application à la SOLVE and INTEG.
In time (already in beta), also for the DM42 ;-)
General Usage Instructions
- Enter a program that defines the differential equation
- Run (XEQ) ODE, select the program
- If applicable, set parameter values, then select the independent variable
- Select the dependent variable
- a. Enter the initial values
b. Enter the value of the independent variable to solve for
c. Set the desired relative accuracy
d. Press the [Yn] key to determine the value of the dependent variable to the desired accuracy.
Solving a first-order ODE
Suppose we want to solve the following ODE
\( \frac {dy}{dx} = e^x + y, x_0=1, y_0=e \)
- Write the function (any alphanumeric name will do)
Code:
00 { 20-Byte Prgm }
01▸LBL "FXY"
02 MVAR "X"
03 MVAR "Y"
04 RCL "X"
05 E^X
06 RCL+ "Y"
07 END
- XEQ "ODE", select the function you just wrote
[FXY]
- Select the independent variable twice..
[X] [X]
(Due to a quirk in the 42S, faithfully reproduced in Free42, the first button selected in a VARMENU will always store whatever value is in X into the corresponding variable; only then can you select a variable without storing anything, resuming execution. That's why the first one has to be pressed twice. Notice that the INTEG menu does not have this, but the SOLVE menu has.)
- Select the dependent variable
[Y]
- a. Enter the initial values for X0 and Y0 :
1 [X0] E^X [Y0]
b. Enter the value of the independent variable you want to solve for :
2 [Xn]
c. Set the desired relative accuracy :
0.00001 [ACC]
d. Press [Yn] to see Yn=14.778109205. The exact value would be 2*e^2 = 14.7781121979
A second example is a hydraulic equation from Gilberto Urroz, probably featured in one of his books that I don't have ;-)
The 48/49-style equation is this:
(1-Q^2/g*((b+2*m*Y)/((b+m*Y)*Y)^3))/(So-(Q*n/Cu)^2*((b+2*Y*sqrt(1+m^2))^(4/3)/((b+m*Y)*Y)^(10/3)))
or
\( { \huge \frac{dX}{dY} = \frac{1 - \frac{Q^2}{g} \cdot \frac{b + 2 \cdot m \cdot Y}{\left( \left( b + m \cdot Y \right) \cdot Y \right) ^3}}{So - \left( \frac{Q \cdot n}{Cu} \right)^2 \cdot \frac { \left( b + 2 \cdot Y \cdot \sqrt{1 + m^2} \right) ^ \frac {4}{3} }{ \left( \left( b + m \cdot Y \right) \cdot Y \right) ^ \frac {10}{3}}} } \)
Mind you, in this case X is the dependent variable and Y the independent
parameter values:
Q: 7
g: 9.81
b: 2
m: 0
So: 0.0003
n: 0.013
Cu: 1
initial values:
Y0 = 2.7
Yn = 3.1
X0 = 0
Xn = ?
- Write the function
Since this time we have included the parameters in the function definition, we can supply their values there. Selecting the independent variable will then only require a single keypress.
This program makes use of the unlimited stack, as it needs 6 stack levels.
Code:
00 { 123-Byte Prgm }
01▸LBL "GU"
02 MVAR "Y"
03 MVAR "X"
04 MVAR "Q"
05 MVAR "g"
06 MVAR "b"
07 MVAR "m"
08 MVAR "So"
09 MVAR "n"
10 MVAR "Cu"
11 FC? 80
12 LNSTK
13 1
14 RCL "Q"
15 X^2
16 RCL÷ "g"
17 RCL "m"
18 RCL× "Y"
19 ENTER
20 RCL+ "b"
21 STO+ ST Y
22 RCL× "Y"
23 LSTO "t"
24 3
25 Y^X
26 ÷
27 ×
28 -
29 RCL "So"
30 RCL "Q"
31 RCL× "n"
32 RCL÷ "Cu"
33 X^2
34 RCL "m"
35 X^2
36 1
37 +
38 SQRT
39 RCL× "Y"
40 STO+ ST X
41 RCL+ "b"
42 0.75
43 1/X
44 Y^X
45 RCL "t"
46 0.3
47 1/X
48 Y^X
49 ÷
50 ×
51 -
52 ÷
53 END
- XEQ "ODE", select the function you just wrote
[GU]
- Set the parameters
7 [Q]
9.81 [g]
2 [b]
0 [m]
0.0003 [So]
0.013 [n]
1 [Cu]
Select the independent variable (only once now since we've entered parameter values..)
[Y]
- Select the dependent variable
[X]
- a. Enter the initial values
2.7 [Y0]
0 [X0]
b. Enter the value of the independent variable you want to solve for
3.1 [Yn]
c. Set the desired relative accuracy
0.0001 [ACC]
d. Press [Xn] to see -8127.82402093
FYI: the 49G took over 5 minutes to come up with an answer (-8127.801629). I can't test it on the DM42 yet, but Free42 on my phone returned the answer immediately.
Solving a system of first-order ODEs
Let's solve the Generalised Lotka-Volterra equation (see the Swissmicros Forum post
here)
- Setup the matrices
"M" 3x1 [[100][ 10][ 1]]
"R" 3x1 [[ 1 ][-1.5][-1.9]]
"A" 3x3 [[ 0 -0.1 0 ][ 0.008 0 -0.1 ] [ 0 0.02 0 ]]
- Write the function
\( \frac {dM}{dt} = DIAG(M) \cdot (R + A \cdot M ) \)
Code:
00 { 32-Byte Prgm }
01▸LBL "GLV"
02 MVAR "t"
03 MVAR "M"
04 RCL "M"
05 XEQ "DIAG"
06 RCL "A"
07 RCL× "M"
08 RCL+ "R"
09 ×
10 END
It is important to note that while t is not used in the definition, it must be declared with MVAR.
- XEQ "ODE", select the function you just wrote
[GLV]
- Select the independent variable twice..
[t] [t]
- Select the dependent variable
[M]
- a. Enter the initial values
0 [t0]
RCL "M" [M0]
b. Enter the value of the independent variable you want to solve for
1 [tn]
c. Set the desired relative accuracy
0.001 [ACC]
d. Press [Mn] to see
Mn = [[133.2678] [5.1880] [0.1725]]
Solving higher-order linear ODEs
As long as the equations are linear, we can rewrite them in vector notation.
Suppose we want to solve the following equation:
x" = -18.75x - 1.962x'
or
\( \frac {d^2 x}{dt^2} = -18.75 \cdot x - 1.962 \cdot \frac {dx}{dt} \)
with initial conditions
t=0 then x=0, x'=6
we can solve this equation rewriting it with a vector:
w' = B*w
with
w = [[ x ][ x' ]], B = [[ 0 1 ][ -18.75 -1.962 ]]
w0 = [[ 0 ][ 6 ]]
- Create the matrices
B 2x2 [[ 0 1 ][ -18.75 -1.962 ]]
w 2x1 [[0][6]]
- Write the function
Code:
00 { 19-Byte Prgm }
01▸LBL "OD2"
02 MVAR "t"
03 MVAR "w"
04 RCL "B"
05 RCLx "w"
06 END
- XEQ "ODE", select the function you just wrote
[OD2]
- Select the independent variable twice
[t] [t]
- Select the dependent variable
[B]
- a. Enter the initial values
RCL "w" [w0]
0 [t0]
b. Enter the value of the independent variable you want to solve for
2 [tn]
c. Set the desired relative accuracy
0.001 [ACC]
d. Press [wn] to see [[1.67168E-1][-6.27E-1]]
Method
- 4th order Runge-Kutta
- Successively halving the step size till the desired relative accuracy is achieved.
for Yn real, test \( ACC < \left| \frac {Y_{n+1} - Y_n}{Y_n} \right| \)
for Yn matrix, test \( ACC < MAX \left( \left| \frac {Y_{n+1}(i) - Y_n(i)}{Y_n(i)} \right| \right) \)
- implemented with a regular MENU i.o. a VARMENU because
- with a VARMENU ODE would also show up in the PGMMENU
- MENU allows defining variables dynamically: X0, t0 etc.
Non-interactive usage, and calling from within a program
Here, the step size and number of steps must be specified, so no automatic adjustment till a desired accuracy is achieved.
- Store the names of the dependent and independent variables in reserved variables ODEP and ODIN, respectively. When they are Y and X you can skip this step, provided ODEP and ODIN do not exist.
- T: n
Z: h
Y: Y0
X: X0
A: ODE function name
XEQ "RK4"
returns
T: n
Z: h
Y: Yn
X: Xn = X0+n*h
A: ODE function name
so you can make a subsequent call to RK4 right away.
the first example would be:
CLV "ODIN"
CLV "ODEP"
A: "FXY"
T: 100
Z: 0.01
Y: 1 E^X
X: 1
XEQ "RK4"
T: 100
Y: 0.01
Z: 14.7781121958
X: 2
XEQ "RK4"
T: 100
Y: 0.01
Z: 60.2566107567
X: 3
the Lotka-Volterra equation:
"t"
ASTO "ODIN"
"M"
ASTO "ODEP"
A: "GLV"
T: 10
Z: 0.1
Y: RCL "M"
X: 0
XEQ "RK4"
T: 10
Z: 0.1
Y: [[133.2681] [5.1880] [1.725E-1]]
X: 1
Code listings
All programs and examples are collected in the attached zip file;
the code listings of the two main programs and the auxiliary one follow here:
The interactive application:
Code:
00 { 353-Byte Prgm }
01▸LBL "ODE"
02▸LBL 12
03 PGMMENU
04 "Select ODE Prog"
05 AVIEW
06 RTN @ will remove local variables
07 CLMENU
08 LASTO ".F"
09 VARMENU IND ".F"
10 "Select indep"
11 ├" Var TWICE"
12 PROMPT
13 LASTO "ODIN"
14 ├"0"
15 LASTO "0X"
16 KEY 1 GTO 01
17 "Select dep. Var"
18 PROMPT
19 LASTO "ODEP"
20 ├"0"
21 LASTO "0Y"
22 KEY 2 GTO 02
23 CLX
24 LSTO ".n"
25 LSTO ".h"
26 CLA
27 ARCL "ODIN"
28 ├"n"
29 LASTO "nX"
30 KEY 3 GTO 03
31 "ACC"
32 KEY 4 GTO 04
33 CLA
34 ARCL "ODEP"
35 ├"n"
36 LASTO "nY"
37 KEY 5 GTO 05
38 KEY 9 GTO 09
39▸LBL 10
40 MENU
41 STOP
42 GTO 10
43▸LBL 01
44 "0X"
45 GTO 00
46▸LBL 04
47 "ACC"
48 GTO 01
49▸LBL 02
50 "0Y"
51 GTO 00
52▸LBL 03
53 "nX"
54▸LBL 00
55 ASTO ST L
56 CLA
57 ARCL IND ST L
58▸LBL 01
59 ASTO ST L
60 SF 25
61 ARCL IND ST L
62 FS?C 25
63 FS? 22
64 LSTO IND ST L
65 CF 22
66 VIEW IND ST L
67 GTO 10
68▸LBL 05
69 CLA
70 ARCL ".F"
71 1
72 STO ".n"
73 RCL IND "nX"
74 RCL IND "0X"
75 -
76 STO ".h"
77 RCL IND "0Y"
78 LASTX
79 XEQ "RK4"
80 X<>Y
81 LSTO IND "nY"
82▸LBL 11
83 2
84 STO÷ ".h"
85 STO× ".n"
86 RCL ".n"
87 RCL ".h"
88 RCL IND "0Y"
89 RCL IND "0X"
90 XEQ "RK4"
91 X<>Y
92 ENTER
93 X<> IND "nY"
94 STO- ST Y
95 MAT?
96 GTO 00
97 X=0?
98 SIGN
99 ÷
100 ABS
101 GTO 01
102▸LBL 00
103 1/X
104 XEQ "DIAG"
105 X<>Y
106 ×
107 RNRM
108▸LBL 01
109 RCL "ACC"
110 X<Y?
111 GTO 11
112 RCL IND "nY"
113 VIEW IND "nY"
114 GTO 10
115▸LBL 09
116 EXITALL
117 GTO 12
118 END
The core, Runge-Kutta 4th order
k0 = h/2.f(xi,yi)
k1 = h/2.f(xi+h/2,yi+k0)
k2 = h/2.f(xi+h/2,yi+k1)
k3 = h/2.f(xi+h,yi+2.k2)
yi+1 = yi + (k0+2.k1+2.k2+k3)/3
To avoid name clashes, numbered registers (in a local REGS) were used when possible
R00 - F (alpha) function name
R01 - X (alpha) indep var name
R02 - Y (alpha) dep var name
R03 - n (integer) number of steps
R04 - h (real) step size
that leaves
"0Y" - initial value dep var - may be a matrix!
"+Y" - deltaY sum - may be a matrix!!
Code:
X Y Z T A
In: x0 y0 h n Function Name
Out: xn yn h n Function Name
00 { 159-Byte Prgm }
01▸LBL "RK4"
02 XEQ 14
03 LSTO "REGS"
04 CLX
05 LASTX
06 FUNC 22
07 ASTO 00
08 SF 25
09 CLA
10 ARCL "ODIN"
11 FC?C 25
12 "X"
13 ASTO 01
14 SF 25
15 CLA
16 ARCL "ODEP"
17 FC?C 25
18 "Y"
19 ASTO 02
20 LSTO IND 01
21 R↓
22 LSTO "0Y"
23 LSTO IND 02
24 R↓
25 2
26 ÷
27 STO 04
28 R↓
29 STO 03
30 LSTO "+Y"
31▸LBL 02
32 CLX
33 XEQ 11
34 STO "+Y"
35 XEQ 10
36 STO+ "+Y"
37 STO+ "+Y"
38 XEQ 11
39 STO+ ST X
40 STO+ "+Y"
41 XEQ 10
42 RCL+ "+Y"
43 3
44 ÷
45 STO+ "0Y"
46 DSE 03
47 GTO 02
48 RCL "0Y"
49 RCL IND 01
50 CLA
51 ARCL 00
52 RTN
53▸LBL 10
54 RCL 04
55 STO+ IND 01
56 X<>Y
57▸LBL 11
58 RCL+ "0Y"
59 STO IND 02
60 XEQ IND 00
61 RCL 04
62 ×
63 RTN
64▸LBL 14
65 FUNC 11
66 5
67 1
68 NEWMAT
69 END
auxiliary routine that turns a 1xn or nx1 vector into a diagonal matrix:
Code:
00 { 37-Byte Prgm }
01▸LBL "DIAG"
02 LSTO "."
03 DIM?
04 +
05 ENTER
06 DSE ST X
07 DIM "."
08 X<> "."
09 TRANS
10 X<> "."
11 STO ST Y
12 DIM "."
13 -
14 X<> "."
15 END
Hope you like it,
Werner
Search
Member List
Calendar
Help
