Casio Classpad fx-CP400: Solving a System of Linear Differential Equations Using Laplace Transforms

Casio Classpad fx-CP400: Solving a System of Linear Differential Equations Using Laplace Transforms

Introduction

The following procedure should work with the Casio Classpad family of calculators (300, 330, fx-CP400, fx-CP500). 

An attempt to use the Casio Classpad to solve the system of differential equations:

dxf/dt = A * xf(t) + B * yf(t) + g1(t)

dyf/dt = C * xf(t) + D * yf(t) + g2(t)

Initial conditions:  xf(0) = xf0,   yf(0) = yf0

The independent variable in all the functions is t.  The Casio Classpad's CAS functions laplace and invLaplace are used to accomplish this task.

Syntax

Laplace Transform of f(t):

laplace(f(t), t, s)

Laplace Transform of a Differential Equation:

laplace(diff eq, independent variable, dependent variable, s)

Inverse Laplace Transformation:

invLaplace(ℒ(s), s, t)

t:  parameter of the original function f(t)

s:  parameter of the transformed function ℒ(s)

The definition of a Laplace Transformation:

ℒ(f(t))(s) = ∫( f(t) * e^(-s*t) dt, t = 0 to ∞)

Now, even though we are working a system of differential equations, the laplace and invLapace commands can only work on function at a time.  Any other variable will be treated as a constant, even it was meant to represent another functions.  The substitution steps will account for this.  

The variables used to respect the Classpad's system variables.  

Here are the steps I used to solve:

dxf/dt = A * xf(t) + B * yf(t) + g1(t)

dyf/dt = C * xf(t) + D * yf(t) + g2(t)

Initial conditions:  xf(0) = xf0,   yf(0) = yf0

Spaces are added to readability.  

1.  Store dxf/dt in the variable f1.   Use xf' to represent dxf/dt.  

2.  Store dyf/dt in the variable f2.   Use yf' to represent dyf/dt.

3.  Execute the following and substitutions:

laplace( f1, t, xf, s ) ⇒ g1

g1 | xf(0) = xf0 and  yf / s = Mp ⇒ g1

4.  Execute the following and substitutions:

laplace( f2, t, yf, s ) ⇒ g2

g2 | yf(0) = yf0 and Lp = Mp ⇒ g2

g2 | xf /s = Lp ⇒ g2

Note Lp = ℒ( xf(t) ) and Mp = ℒ( yf(t) )

5.  Now solve the system of transformed equations: 

solve( {g1, g2}, {Lp, Mp} ) ⇒ lists

6.  Take the inverse laplace transforms:

invLaplace( getRight( lists[1], s, t ) ) ⇒ h1

invLaplace( getRight( lists[2], s, t ) ) ⇒ h2

The solutions are:

xf(t) = h1

yf(t) = h2

The getRight command extracts the right side of an equation.

Examples

Example 1:

xf ' = -0.08 * xf + 0.02 * yf + 6

yf ' = 0.08 * xf - 0.08 * yf

Initial conditions:  xf(0) = 0, yf(0) = 150

Example 2:

xf ' = 2 * xf + yf + t^2

yf ' = -2 * xf + 5

Initial conditions:  xf(0) = 2, yf(0) = 5

Casio Classpad Program:  laplaced

Download here:  https://drive.google.com/file/d/1zxTEynCjXmBrPU9pYqmXo3kF_so7fiI-/view?usp=share_link

Code:

' setup

SetRadian

' local variables

Local str1, str2, f1, f2

Local g1, g2, lists, h1, h2

Local Mp, Lp, t, s, xval, yval

' main

ClrText

InputStr str1, "xf'(xf,yf,g1(t))="

InputStr str2, "yf'(xf,yf,g2(t))="

Input xval, "xf(0)? "

Input yval, "yf(0)? "

StrJoin "xf'=", str1, str1

StrJoin "yf'=", str2, str2

strToExp(str1) ⇒ f1

strToExp(str2) ⇒ f2

' transform and solve

laplace(f1, t, xf, s) ⇒ g1

g1 | xf(0) = xval and yf/s = Mp ⇒ g1

laplace(f2, t, yf, s) ⇒ g2

g2 | yf(0) = yval and Lp = Mp ⇒ g2

g2 | xf/s = Lp ⇒ g2

Print "Laplace Transforms:"

PrintNatural g1

Print g1

PrintNatural g2

Print g2

solve({g1, g2},{Lp, Mp})⇒lists

invLaplace(getRight(lists[1]),s,t) ⇒ h1

invLaplace(getRight(lists[2],s,t) ⇒ h2

Print "Solutions: "

Print "xf = "

PrintNatural h1

Print h1

Print "yf = "

PrintNatural h2

Print h2

Note:  PrintNatural displays the expression in textbook form in a popup box, but it does not add anything to the text output terminal.  This is why I included both PrintNatural and Print commands.  

Source

Kreyszig, Erwin.  Advanced Engineering Mathematics  8th Edition  John Wiley & Sons, Inc:  New York, NY.  1999.  ISBN 0-471-15496-2

All original content copyright, © 2011-2023.  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-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.

Differential Equations #12: Reduction of Order

Introduction

The general n-order differential equation

can be translated into a system of first-order differential equations:

For a second-order equation

y" + A*y' + B*y = f(x)

The system is

y'_1 = y_2
y'_2 = f(x) - A*y_2 - B*y_1

Now time for examples.

1. y" + 3*y' - 4*y = 2

Let y_1 = y and y_2 = y'_1. Then:
y_2 = y', and

y'_2 = y"
y'_2 = 2 - 3*y' - 4*y
= 2 - 3*y_2 - 4*y_1

So our system is:
y'_1 = y_2
y'_2 = 2 - 3 * y_2 - 4 * y_1


2. 3*y" + 6*y' = -9

Let y_1 = y and y_2 = y'. Then:

y'_1 = y_2
y'_2 = y"
3*y" = -9 - 6*y'
y" = -3 - 2*y'
y'_2 = -3 - 2*y_2

The system is:
y'_1 = y_2
y'_2 = -3 - 2*y_2


We can go the other way, from a system of linear differential equations to a single multi-order differential equation.

3. y'_1 = y_2
y'_2 = y_1 + y_2

Again, let y_1 = y and y_2 = y'. Then:
y'_2 = y_1 + y_2
y" = y + y'
y" - y' - y = 0

4. y'_1 = y_2
y'_2 = 2/7 * y_1 - 3/7 * y_2 + 1/7

With y_1 = y and y_2 = y',
y'_2 = 2/7 * y_1 - 3/7 * y_2 + 1/7
y" = 2/7 * y - 3/7 * y' + 1/7
y" + 3/7 * y' - 2/7 * y = 1/7
7*y" + 3*y' - 2*y = 1

This concludes this blog post of Differential Equations. In a week I will be headed to Fort Collins. CO for the HHC 2013 conference. I am so excited - may the bad weather go away and everyone be safe!

Eddie


This blog is property of Edward Shore. 2013

Differential Equations #10: Critical Points

Critical Points and Determining What Happens

In this blog entry we are working with a system of two equations:

x' = f(x,y)
y' = g(x,y)

where x and y are functions of a independent variable, say t for example. Well treat t as a time variable.

Today's blog will cover a three step process:
1. Finding Critical Points
2. Determining the Jacobian Matrix
3. Finding the Eiegenvalues of that Jacobian Matrix

1. Finding Critical Points

It is relatively easy to find the critical points of a system. Add that needs to be done is to set x' = 0 and y' = 0. The next step is to solve for x and y.

2. Determining the Jacobian Matrix

The Jacobian matrix of the system is

---

Hint: If the system is
x' = Ax + By
y' = Cx + Dy

Then the system's critical point is (0, 0) and the Jacobian Matrix is

J = [ [A, B] [C, D] ]

---

3. Determine the Eigenvalues

As we recall the eigenvalues of a 2 × 2 matrix can be determined by solving the following characteristic equation:

λ^2 - trace * λ + determinant = 0

Where:
* trace is the sum of the matrix's diagonal elements
* determinant is the determinant of the matrix

The roots are crucial and determine the behavior of the critical point.

Types of critical points:
* source: individual curve solutions x(t) and y(t) are trajectories going away from the critical point
* sink: individual curve solutions x(t) and y(t) are trajectories go towards the critical point
* saddle point: the critical point acts as a sink for some trajectories and a source for other trajectories
* center: trajectories orbit around the critical point, most likely circular or elliptical orbits

Determining the type of critical points:
* λ1 and λ2 are real and positive: critical point is a source.
* λ1 and λ2 are real and negative: critical point is a sink
* λ1 and λ2 are real and have opposite signs: critical point is a saddle
* λ is a double root and it's positive: critical point is a source
* λ is a double root and it's negative: critical point is a sink
* λ = S ± Ti and S is positive: point is a source, trajectories are spiral
* λ = S ± Ti and S is negative: point is a sink, trajectories are spiral
* λ = ± Ti: critical point is a center

Source: http://kaharris.org/teaching/216/Lectures/lec28/lec28.pdf

Let's work with some examples.

1.
x' = x + 4y
y' = 2x + 3y


Critical Point:
0 = x + 4y
0 = 2x + 3y

x = -4y
0 = 2(-4y) + 3y which implies y = 0 and x = 0

Our critical point is (0,0). The Jacobian matrix is J = [ [1, 4] [2, 3] ].

trace(J) = 1 + 3 = 4
det(J) = 1 * 3 - 2 * 4 = -5

The eigenvalues are:
λ^2 - 4*λ - 5 = 0
(λ - 5)(λ + 1) = 0
λ1 = 5, λ2 = -1
Real and opposite signs, the critical point (0,0) is a saddle point.

2.
x' = x + 2y - 6
y' = 6x - 3y + 24

Critical Point:
0 = x + 2y - 6
0 = 6x - 3y + 24

Which leaves is with the system of
x + 2y = 6
6x - 3y = -24

Solving the system yields the critical point (-2, 4).

With
f = x + 2y - 6
df/dx = 1
df/dy = 2
g = 6x - 3y + 24
dg/dx = 6
dg/dy = -3,

The Jacobian matrix is
J = [ [1, 2] [6, -3] ]

The eigenvalues of J are 3 and -5. Since the eigenvalues are real and are of opposite signs, the critical point (-2, 4) is a saddle point.

3.
x' = 2x + 3y
y' = x - 2y^2

Finding the critical points:
0 = 2x + 3y
0 = x - 2y^2

x = 2y^2
0 = 2(2y^2) + 3y
0 = 4y^2 + 3y
0 = y * (4y + 3)
which implies that we have two critical points, one with y = 0 and the other with y = -3/4.

If y = 0, then x = 2*0^2 = 0. One critical point is (0,0).
If y = -3/4, then x = 2*(-3/4)^2 = 9/8. The other critical point is (9/8, -3/4).

With
f = 2x + 3y, df/dx = 2, df/dy = 3
g = x - 2y^2, dg/dx = 1, dg/dy = -4y

Critical Point (0,0):
J = [ [2, 3] [1, 0] ]
And the eigenvalues of J are -1 and 3.
The critical point (0,0) is a saddle point.

Critical Point (9/8, -3/4):
J = [ [2, 3] [1, 3] ]
And the eigenvalues are (5 ± √13)/2. Approximate values are 4.302776 and 0.697224. The critical point (9/8, -3/4) is a source.

Hope this helps and we'll talk next time. Please feel free to post questions and comments.

Eddie


This blog is property of Edward Shore. 2013