(Free42) interactive ODE
03-31-2021, 06:39 PM (This post was last modified: 04-09-2021 07:35 AM by Werner.)
Post: #1
 Werner Senior Member Posts: 606 Joined: Dec 2013
(Free42) interactive ODE
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

1. Enter a program that defines the differential equation
2. Run (XEQ) ODE, select the program
3. If applicable, set parameter values, then select the independent variable
4. Select the dependent variable
5. 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$$

1. 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
2. XEQ "ODE", select the function you just wrote
[FXY]
3. 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.)
4. Select the dependent variable
[Y]
5. 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 = ?

1. 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
2. XEQ "ODE", select the function you just wrote
[GU]
3. 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]
4. Select the dependent variable
[X]
5. 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)

1. 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 ]]
2. 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.
3. XEQ "ODE", select the function you just wrote
[GLV]
4. Select the independent variable twice..
[t] [t]
5. Select the dependent variable
[M]
6. 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 ]]

1. Create the matrices
B 2x2 [[ 0 1 ][ -18.75 -1.962 ]]
w 2x1 [[0][6]]
2. 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
3. XEQ "ODE", select the function you just wrote
[OD2]
4. Select the independent variable twice
[t] [t]
5. Select the dependent variable
[B]
6. 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)$$
- 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.

1. 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.
2. 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

Attached File(s)