HP Prime and TI-84 Plus: Runge Kutta 4th Order

HP Prime and TI-84 Plus:  Runge Kutta 4^th Order

The Runge Kutta 4^th Order is a method for solving differential equations involving the form:  dy/dx = f(x,y), where:

x_n+1 = x_n + h
y_n+1 = y_n + (k1 + 2*k2 + 2*k3 + k4)/6 

Where:

k1 = h * f(x_n, y_n)
k2 = h * f(x_n + h/2, y_n + k1/2)
k3 = h * f(x_n + h/2, y_n + k2/2)
k4 = h * f(x_n + h, y_n + k3)

Variables used:

A = x_n
B = y_n
C = x_n+1
D = y_n+1
H = step
K = k1
L = k2
M = k3
N = k4

Results are stored in lists L1 and L2.  L1 represents the x coordinates, L2 represents the y coordinates.

This adopts the RK4 method that I programmed for the Casio fx-5800p, which you can see here:  http://edspi31415.blogspot.com/2015/01/fx-5800p-runge-kutta-4th-order-method.html

The screen shots show the example dy/dx = x^2 + y^2/3, with initial conditions y(0) = 1.  10 steps are plotted with step size 0.1.  

HP Prime:  DIFFTBL and RK4

Both DIFFTBL and RK4 return the same output.  The difference is that DIFFTBL uses an Input box and no-pass through arguments and RK4 uses five pass-through arguments.  This way, RK4 could be used as a subroutine.

Both return the results in a matrix, M1.  Rev. 7820 firmware is used.

DIFFTBL:

EXPORT DIFFTBL()

BEGIN

// EWS 2015-05-29

// Radian

HAngle:=0;

// Input

LOCAL f;

  

INPUT({{f,[8]},A,B,H,S},"Data Input",

{"dy/dx=","x0=","y0=","Step=",

"# Steps="});

L1:={A}; L2:={B};

// Main Loop

FOR I FROM 1 TO S DO

// k1

X:=A;

Y:=B;

K:=H*EVAL(f);

// k2

X:=A+H/2;

Y:=B+K/2;

L:=H*EVAL(f);

// k3

Y:=B+L/2;

M:=H*EVAL(f);

// k4

X:=A+H;

Y:=B+M;

N:=H*EVAL(f);

// point

A:=A+H;

B:=B+(K+2*L+2*M+N)/6;

L1:=CONCAT(L1,{A});

L2:=CONCAT(L2,{B});

END;

M1:=list2mat(CONCAT(L1,L2),

SIZE(L1));

M1:=TRN(M1);

RETURN M1;

END;

RK4:

EXPORT RK4(f,A,B,H,S)

BEGIN

// EWS 2015-05-29

// Radian

HAngle:=0;

// Input

L1:={A}; L2:={B};

// Main Loop

FOR I FROM 1 TO S DO

// k1

X:=A;

Y:=B;

K:=H*EVAL(f);

// k2

X:=A+H/2;

Y:=B+K/2;

L:=H*EVAL(f);

// k3

Y:=B+L/2;

M:=H*EVAL(f);

// k4

X:=A+H;

Y:=B+M;

N:=H*EVAL(f);

// point

A:=A+H;

B:=B+(K+2*L+2*M+N)/6;

L1:=CONCAT(L1,{A});

L2:=CONCAT(L2,{B});

END;

M1:=list2mat(CONCAT(L1,L2),

SIZE(L1));

M1:=TRN(M1);

RETURN M1;

END;

TI-84 Plus:  DIFFPLOT

DIFFPLOT plots the results in a statistics xy-plot.  The equation variable Y0 is used, but subsequently deleted.  I left out any color commands to make this easily adoptable for both the monochrome and color TI-84 Plus calculators. 

Do not forget to enclose dy/dx in quotes, unless an input error occurs

Radian:FnOff :PlotsOff

Disp "DIFF. EQUATION"

Disp "Y₀=DY/DX"

Disp "USE QUOTES"

Input Y₀

Input "X0=",A

Input "Y0=",B

Input "STEP=",H

Input "NO. STEPS=",S

{A}→L₁

{B}→L₂

For(I,1,S)

A→X:B→Y

H*Y₀→K

A+H/2→X:B+K/2→Y

H*Y₀→L

B+L/2→Y

H*Y₀→M

A+H→X:B+M→Y

H*Y₀→N

A+H→A

B+(K+2*L+2*M+N)/6→B

augment(L₁,{A})→L₁

augment(L₂,{B})→L₂

End

FnOff 0

PlotsOn 1

Plot1(xyLine,L₁,L₂)

ZoomStat

This blog is property of Edward Shore.  2015.

fx-5800p: Runge-Kutta 4th Order Method - Differential Equations

Runge-Kutta 4th Order Method - Differential Equations

The program RK4 estimates solutions to the differential equation y dy/dx = f(x,y) given the initial condition (x_n, y_n). With step size h, the point (x_n+1, y_n+1) is calculated by:

x_n+1 = x_n + h
y_n+1 = y_n + (k1 + 2*k2 + 2*k3 + k4)/6

Where:

k1 = h * f(x_n, y_n)
k2 = h * f(x_n + h/2, y_n + k1/2)
k3 = h * f(x_n + h/2, y_n + k2/2)
k4 = h * f(x_n + h, y_n + k3)

The program RK4 uses a subroutine FXY. The program FXY is where the y'(x) is stored in terms of X and Y. The result is stored in the variable Z.

Variables used in RK4 and FXY:
A = x_n
B = y_n
C = x_n+1
D = y_n+1
H = step
K = k1
L = k2
M = k3
N = k4
E, X, Y, and Z are also used

Remember for the fx-5800p, all variables (A-Z, and Z[#]s) are global, meaning they will carry over from programs and subroutines.

Anything following the double slash (//) is a comment.

Once the first point (x_n+1, y_n+1) is calculated, RK4 will ask if you want the next point (x_n+2, y_n+2) calculated. Enter 1 at the prompt for "Yes".

Program FXY
type f(X,Y) here → Z

// A RETURN command is implicit.

Program RX4
"INITIAL COND."
"X0"?→ A
"Y0"?→ B
"STEP"?→ H
Lbl 0 // main routine
A → X // k1
B → Y
Prog "FXY"
H * Z → K
A + H ÷ 2 → X // k2
B + K ÷ 2 → Y
Prog "FXY"
H * Z → L
B + L ÷ 2 → Y // k3
Prog "FXY"
H * Z → M
A + H → X // k4
B + M → Y
Prog "FXY"
H * Z → N
X → C // x_n+1
B + (K + 2*L + 2*M + N) ÷ 6 → D // y_n+1
"RESULT"
"X1 = "
C ◢
"Y1 = "
D ◢
"NEXT POINT?"
"1 = YES 0 = NO"
? → E
E ≠ 1 ⇒ Goto E
C → A
D → B
Goto 0
Lbl E
"DONE"



Example: y'(x) = y - x with the initial condition (0,2).
Find y when x = 0.1 and x = 0.2, respectively.
In this case, our step is h = 0.1.

Program FXY will have this:
Y - X → Z


Running RK4 gives these results:
x = 0.1, y = 2.205170833
x = 0.2, y = 2.421402571

Eddie

This blog is property of Edward Shore. 2015

Differential Equations #11: Numerical Solutions with the Runge-Kutta 4th Order Method

Update:  One of the numerical example (example 1) has incorrect numerical results.  This blog has been edited.  The results posted are now correct. Apologizes. -  Eddie

1/4/2015

Introduction

In this blog we will explore one of the ways to solve the differential equation y' = f(x,y) by numerical methods. One of the common methods is the Runge-Kutta 4th Order Method, probably the most common methods.

With the initial condition (x_n, y_n) and step size h, the next point (x_n+1, y_n+1) is calculated by the following equations:
x_n+1 = x_n + h
y_n+1 = y_n + (k1 + 2*k2 + 2*k3 + k4)/6 with

k1 = h * f(x_n, y_n)
k2 = h * f(x_n + h/2, y_n + k1/2)
k3 = h * f(x_n + h/2, y_n + k2/2)
k4 = h * f(x_n + h, y_n + k3)

(Almost) obviously use either a calculator or mathematical software for this process.

Examples:
1. dy/dx = -0.23 * (y + 1)
Initial Condition: (0, 24)
We want y1 when x1 = 1 and y2 when x2 = 2.

Observe that h = 1 with x0 = 0 and y0 = 24. Note the lack of x term in f(x,y).

Then: x1 = 0 + 1 = 1
k1 ≈ -5.75
k2 ≈ -5.08875
k3 ≈ -5.16479375
k4 ≈ -4.562097438
y1 ≈ 18.86346918

x2 = 1 + 1 = 2
k1 ≈ -4.568597911
k2 ≈ -4.043209151
k3 ≈ -4.103628858
k4 ≈ -3.624763273
y2 ≈ 14.78229631

So our points are:
(0, 24)
(1, 18.86346918)
(2, 14.78229631)

2. y' = y - x with I.C. y(0) = 2
Find y(0.1) and y(0.2).

Note h = 0.1 with x0 = 0, y0 = 2, and f(x,y) = -x + y.

x1 = 0.1
k1 = 0.1 * f(0, 2) = 0.2
k2 = 0.1 * f(0.5 * 0.1, 2 + 0.5 * 0.2) = 0.205
k3 = 0.1 * f(0.5 * 0.1, 2 + 0.5 * 0.205) = 0.20525
k4 = 0.1 * f(0.1, 2 + 0.20525) = 0.210525
y1 = 2 + (0.2 + 2 * 0.205 + 2 * 0.20525 + 0.210525)/6 ≈ 2.205171

x2 = 0.2
k1 ≈ 0.210517
k2 ≈ 0.216043
k3 ≈ 0.216319
k4 ≈ 0.222149
y2 ≈ 2.421203

So our points are:
(0, 2)
(0.1, 2.205171)
(0.2, 2.421203)


Talk to you all next time!
Eddie

This blog is property of Edward Shore. 2013

3rd Order Runge-Kutta - HP 15C

This program uses a 3rd Order Runge-Kutta method to assist in solving a first order-differential equation.

Given the initial condition (x0, y0) to the differential equation:

dy/dx = f(x, y)

Find the value of y1 where h is a given step size (preferably small). The value of y1 is given by the following equations:

y1 = y0 + k1/4 + (3 * k3)/4

Where:
k1 = h * f(x0, y0)
k2 = h * f(x0 + h/3, y0 + k1/3)
k3 = h * f(x0 + (2 * h)/3, y0 + (2 * k2)/3)

Error estimated on the order of h^4

Source: Smith, Jon M. Scientific Analysis on the Pocket Calculator. John Wiley & Sons, Inc.: New York 1975 pg 174

Memory Registers Used

R0 = h
R1 = x
R2 = y
R3 = x1 (result)
R4 = y1 (result)
R5 = x0 (copied from R1)
R6 = y0 (copied from R2)
R7 = k1
R8 = k2
R9 = k3

Labels Used

Label A = Main Program
Label 0 = Program where equation is maintained

Instructions

1. Store h, x0, and y0 in memory registers R0, R1, and R2 respectively.
2. Enter or edit the equation in program mode using Label 0. Use RCL 1 for x and RCL 2 for y.
3. Run Program A. The first answer is x1. Press [R/S] for y1.
4. To find (x2, y2), store x1 in R1 and y1 in R2 and run Program A again. You can do this by pressing [STO] [ 2 ] [x<>y] [STO] [ 1 ] [ f ] [ √ ] (A) immediately following program execution.

There are 59 program steps in the main program.

Key		Code
LBL A		42	21	11
RCL 1			45	1
STO 5			44	5
RCL+ 0		45	40	0
STO 3			44	3
RCL 2			45	2
STO 6			44	6
GSB 0			32	0
RCLx 0		45	20	0
STO 7			44	7
RCL 5			45	5
RCL 0			45	0
3				3
÷				10
+				40
STO 1			44	1
RCL 6			45	6
RCL 7			45 	7
3				3
÷				10
+				40
STO 2			44	2
GSB 0			32	0
RCLx 0		45	44	0
STO 8			44	8
RCL 5			45	5
RCL 0			45	0
2				2
x				20
3				3
÷				10
+				40
STO 1			44	1
RCL 6			45	6
RCL 8			45	8
2				2
x				20
3				3
÷				10
+				40
STO 2			44	2
GSB 0			32	0
RCLx 0		45	20	0
STO 9			44	9
RCL 3			45	3
R/S				31
RCL 6			45	6
RCL 7			45	7
4				4
÷				10
+				40
RCL 9			45	9
3				3
x				20
4				4
÷				10
+				40
STO 4			44	4
RTN			43	32

Example 1

dy/dx = x^2 + sin(x * y)

Initial Conditions: (1, 1), let h = 0.01

Program for equation:

 
LBL 0
RAD
RCL 1
x^2
RCL 1
RCLx 2
SIN
+
RTN

Results:

x	y
1.0000	1.0000
1.0100	1.0186
1.0200	1.0375
1.0300	1.0568

Example 2

dy/dx = x^3 * y - y

Initial Condition: (1,1); Step size h = 0.01

Program for the equation:

 
LBL 0
RCL 1
3
y^x
RCLx 2
RCL- 2
RTN

Results:

x	y
1.0000	1.0000
1.0100	1.0002
1.0200	1.0006
1.0300	1.0014

This blog is property of Edward Shore. © 2012