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.

Differential Equations #3: General First-Order Linear Differential Equations

On this installment of the Differential Equation series, we will look at general linear differential equations which takes the form

y' + p(x) * y = q(x)

On this blog entry, we are going to go over a general procedure on how to solve these type of equations. The procedure requires a multiplicative factor called the integrating factor.

Procedure to solve y' + p(x) * y = q(x)

1. Calculate the integrating factor I = e^( ∫ p(x) dx).

2. Multiply the factor, I, to the equation, resulting in:
I * y' + (p(x) * I) * y = q(x) * I

Since I' = d/dx ( ∫ p(x) dx) * e^( ∫ p(x) dx) = p(x) * e^( ∫ p(x) dx), I' = I * p(x).

Then I * y' + I' * y = q(x) * I

By the derivative product rule, d/dx ( I * y ) = q(x) * I

3. Take the integral with respect to x on both sides:

I * y = ∫ (q(x) * I) dx + C

4. Solve for y. If an initial condition is given, you can solve for C.

y = ( ∫ q(x) * I dx + C) / I

---

We can summarize this produce by the following:

To solve y' + p(x) * y = q(x)

y = ( ∫ q(x) * I dx + C) / I, where I = e^( ∫ p(x) dx)



Let's work on some examples.

1. y' + 3*y = x

Here p(x) = 3 and q(x) = x. The integrating factor is I = e^( ∫ p(x) dx), so for this example, I = e^( ∫ 3 dx) = e^(3*x). Now:

e^(3*x) * y' + 3 * y * e^(3*x) = x * e^(3*x)
d/dx (y * e^(3*x) ) = x * e^(3*x)
∫ d/dx (y * e^(3*x)) dx = ∫ x * e^(3*x) dx

Using integration by parts on the right side:
y * e^(3*x) = (x*e^(3*x))/3 - (e^(3*x))/9 + C

y = x/3 - 1/9 + C*e^(-3*x)

2. y' + 2/x * y = x^3

Here p(x) = 2/x and q(x) = x^3, and the integrating factor is
I = e^( ∫ 2/x dx) = e^(2*ln x) = e^(ln (x^2)) = x^2

Then:
d/dx (y * x^2) = x^2 * x^3
d/dx (y * x^2) = x^5
∫ d/dx (y * x^2) dx = ∫ x^5 dx
y * x^2 = x^6/6 + C

y = x^4/4 + C*x^(-2)

3. y' + x*y = 2*x, y(0)=3

First, solve for y:
p(x) = x
q(x) = 2*x
I = e^( ∫ x dx) = e^(x^2/2)

d/dx (e^(x^2/2) * y) = e^(x^2/2) * 2 * x
e^(x^2/2) * y = ∫ e^(x^2/2) * 2 * x dx + C
e^(x^2/2) * y = 2*e^(x^2/2) + C
y = 2 + C*e^(-x^2/2)

Using the initial condition y(0)=3 to solve for C:
3 = 2 + C*e^0
C = 1

The final answer is:
y = 2 + e^(-x^2/2)

The next time we will look at (simple) second-order differential equations, and how solving a quadratic equation helps solve the differential equations. Please leave comments and questions below.

Take care, Eddie



This blog is property of Edward Shore. 2013