Saturday, August 1, 2020

Casio fx-9750GIII and fx-CG50: System of Two Differential Equations, Runge Kutta 4th Order

Casio fx-9750GIII and fx-CG50: System of Two Differential Equations, Runge Kutta 4th Order

Introduction




The program TWORK4 uses the Runge Kutta 4th Order method to solve the following system of differential equations:

dx/dt = f(t, x, y)
dy/dt = g(t, x, y)

with initial conditions x0 = x(t0) and y0 = y(t0)

The next step is calculated with step h from:

x1 ≈ x0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6
y1 ≈ y0 + (l1 + 2 * l2 + 2 * l3 + l4) / 6

k1 = h * f(t0, x0, y0)
l1 = h * g(t0, x0, y0)

k2 = h * f(t0 + h / 2, x0 + k1 / 2, y0 + l1 / 2)
l2 = h * g(t0 + h / 2, x0 + k1 / 2, y0 + l1 / 2)

k3 = h * f(t0 + h / 2, x0 + k2 / 2, y0 + l2 / 2)
l3 = h * g(t0 + h / 2, x0 + k2 / 2, y0 + l2 / 2)

k4 = h * f(t0 + h, x0 + k3, y0 + l3)
l4 = h * g(t0 + h, x0 + k3, y0 + l3)

For the next step set t0 = t0 + h, x0 = x1, and y0 = y1

Inputs:

DX/DT:  Enter dx/dt as a function of T, X, and Y.   T is the independent variable.

DY/DT:  Enter dy/dt as a function of T, X, and Y.   T is the independent variable.

T0, X0, Y0:  Enter the initial conditions

STEP:  Enter the step size

ITERATIONS: Enter the number of iterations desired.  This allows you to calculate a far point in one leap.  Example, if your initial condition is to = 0 and you want to find the point when t = 1 using the step size of 0.1, enter 0.1 for STEP and 10 for ITERATIONS.

Casio fx-9750GIII and fx-CG 50 Program: TWORK4

Notes: 

*  Although the code for the two calculators are the same, the programs will need to be programmed on each calculator separately.

*  The slash character (/) is accessed from the CHAR submenu.  The submenu shows up at the top-level menu which would read:
TOP/BOTTOM/SEARCH/MENU/A ←→ a/CHAR

* fn1 and fn2 are stored as function memories. 

"2020-06-18 EWS"
Rad
"DX/DT="? → fn1
"DY/DT="? → fn2
"T0"? → U
"X0"? → A
"Y0"? → B
"STEP?" → H
4 → Dim List 25
4 → Dim List 26
Lbl 0
"ITERATIONS"? → N
For 1 → I To N
A → X: B → Y: U → T
H * fn1 → List 25[1]
H * fn2 → List 26[1]
A + List 25[1] ÷ 2 → X
B + List 26[1] ÷ 2 → Y
U + H ÷ 2 → T
H * fn1 → List 25[2]
H * fn2 → List 26[2]
A + List 25[2] ÷ 2 → X
B + List 26[2] ÷ 2 → Y
H * fn1 → List 25[3]
H * fn2 → List 26[3]
A + List 25[3] → X
B + List 26[3] → Y
U + H → T
H * fn1 → List 25[4]
H * fn2 → List 26[4]
A + (List 25[2] + List 25[3] + Sum List 25) ÷ 6 → A
B + (List 26[2] + List 26[3] + Sum List 26) ÷ 6 → B
U + H → U
Next
ClrText
"(T, X, Y)"
{U, A, B} ◢
Menu "NEXT?", "YES", 0, "NO", 1
Lbl 1
"DONE"

Example

dx/dt = x sin y

dy/dt = y^2/500 + y/3 - y

Initial conditions: x(0) = 0.1, y(0) = 0.1, h = 0.1

T0 = 0
STEP = 0.1

ITERATIONS: 10
Results:  (T, X, Y): (1, 0.1075645239, 0.05134921355)

ITERATIONS: 10
Results:  (T, X, Y): (2, 0.1116714419, 0.02636554462)

Source:

John W. Harris and Horst Stocker.  Handbook of Mathematics and Computational Science.  Springer: New York.  2006.  ISBN 978-0-387-94746-4

Stack Exchange.  "Help with using the Runge-Kutta 4th order method on a system of 2 first order ODE'S"  Asked March 2014.   Last updated February 2019.  https://math.stackexchange.com/questions/721076/help-with-using-the-runge-kutta-4th-order-method-on-a-system-of-2-first-order-od  Retrieved June 17, 2020

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Casio fx-CG50 and Swiss Micros DM32: HP 16C’s Bit Summation

  Casio fx-CG50 and Swiss Micros DM32: HP 16C’s Bit Summation The HP 16C’s #B Function The #B function is the HP 16C’s number of...