TI-84 Plus CE and Numworks: System of Linear Differential Equations (Simple System)
Solving The System
The program presented today will solve the following system of equations:
dx/dt = A * x(t) + B * y(t)
dy/dt = C * x(t) + D * y(t)
with initial conditions x(0) = x0 and y(0) = y0
The formulas used in the program are derived using the Laplace Transform.
dx/dt = A * x(t) + B * y(t)
dy/dt = C * x(t) + D * y(t)
Laplace Transforms to be used:
ℒ{f(t)} = F(s)
ℒ{df/dt} = s * F(s) - f(0)
Apply the Laplace Transform transform:
s * X(s) - x0 = A * X(s) + B * Y(s)
s * Y(s) - y0 = C * X (s) + B * Y(s)
(s - A) * X(s) + (-B) * Y(s) = x0
(-C) * X(s) + (s - D) * Y(s) = y0
Define three matrices:
Md = [ [ s - A, -B ], [ -C, s - D ] ]
where
det(Md) = s^2 - (A + D) * s + (A * D - B * C)
which can be factored into:
(s - P) * (s - Q)
where:
P = (( A + D ) + √(( A + D)^2 - 4 * (A * D - B * C)) /2
Q = (( A + D ) - √(( A + D)^2 - 4 * (A * D - B * C)) /2
Mx = [ [ x0, -B ], [ y0, s - D ] ]
where
det(Mx) = x0 * (s - D) + B * y0
My = [ [ s - A, x0 ], [ - C, y0 ] ]
where
det(My) = (s - A) * y0 + C * x0
Per Cramer's Rule, the solution to X(s) and Y(s) are:
X(s) = det(Mx) / det(Md)
Y(s) = det(My) / det(Md)
What happens next depends of the roots P and Q. We're assuming that P and Q are real here but the situation applies if P and Q are complex.
If P ≠ Q:
If P ≠ Q, then after simplifying and finding partial fractions:
X(s) = ( x0 * (s - D) + B * y0 ) / ((s - P) * (s - Q)) = α / (s - P) + β / (s - Q)
where
α = ( B * y0 + (P - D) *x0) / (P - Q)
β = x0 - α
Y(s) = ( (s - A) * y0 + C * x0 ) / ((s - P) * (s - Q)) = γ / (s - P) + δ / (s - Q)
where
γ =(C * x0 + y0 * (P - A)) / (P - Q)
δ = y0 - γ
Applying the inverse Laplace transform:
x(t) = α * e^(P * t) + β * e^(Q * t)
y(t) = γ * e^(P * t) + δ * e^(Q * t)
If P = Q:
If P = Q, then after simplifying and finding partial fractions:
X(s) = ( x0 * (s - D) + B * y0 ) / (s - P)^2 = α / (s - P)^2 + β / (s - Q)
where
α = B * y0 + (P - D) * x0
β = x0
Y(s) = ( (s - A) * y0 + C * x0 ) / (s - P)^2 = γ / (s - P)^2 + δ / (s - Q)
where
γ =C * x0 + y0 * (P - A)
δ = y0
Applying the inverse Laplace transform:
x(t) = α * t * e^(P * t) + β * e^(Q * t)
y(t) = γ * t * e^(P * t) + δ * e^(Q * t)
TI-84 Plus CE Program: LINSYS1
Download the program here: https://drive.google.com/file/d/14zYpT4xhLRGVvBNriOacIOGBpqXLWswR/view?usp=share_link
Program Listing:
Radian
a+bi
ClrHome
Disp "DX/DT=AX+BY DY/DT=CX+DY"
Prompt A,B,C,D
Input "X0? ",M
Input "Y0? ",N
((A+D)+√((A+D)²-4*(A*D-B*C)))/2→P
((A+D)-√((A+D)²-4*(A*D-B*C)))/2→P
If P≠Q
Then
(B*N-D*M+P*M)/(P-Q)→R
M-R→S
(C*M-A*N+P*N)/(P-Q)→U
N-U→V
ClrHome
Output(2,1,"X=Re^(PT)+Se^(QT)")
Output(3,1,"R= "+toString(R))
Output(4,1,"P= "+toString(P))
Output(5,1,"S= "+toString(S))
Output(6,1,"Q= "+toString(Q))
Pause
ClrHome
Output(2,1,"Y=Ue^(PT)+Ve^(QT)")
Output(3,1,"U= "+toString(U))
Output(4,1,"P= "+toString(P))
Output(5,1,"V= "+toString(V))
Output(6,1,"Q= "+toString(Q))
Pause
Else
B*N-D*M+M*P→P
M→S
C*M-A*N+N*P→U
N→V
ClrHome
Output(2,1,"X=RTe^(PT)+Se^(PT)")
Output(3,1,"R= "+toString(R))
Output(4,1,"P= "+toString(P))
Output(5,1,"S= "+toString(S))
Pause
ClrHome
Output(2,1,"Y=UTe^(PT)+Ve^(PT)")
Output(3,1,"U= "+toString(U))
Output(4,1,"P= "+toString(P))
Output(5,1,"V= "+toString(V))
Pause
End
ClrHome
Python: dfsys1.py
Link from my Numworks page: https://my.numworks.com/python/ews31415/dfsys1
Modules:
cmath: complex math
mathplotlib.pyplot: screen plotting and text placement
Script listing:
from cmath import *
from matplotlib.pyplot import *
# 2023-01-28 EWS
# parameter input
print("systems of differential\nequations")
print("dx/dt=a*x+b*y")
print("dy/dt=c*x+d*y")
a=float(input("a? "))
b=float(input("b? "))
c=float(input("c? "))
d=float(input("d? "))
m=float(input("x0? "))
n=float(input("y0? "))
# characteristic roots
p=((a+d)+sqrt((a+d)**2-4*(a*d-b*c)))/2
q=((a+d)-sqrt((a+d)**2-4*(a*d-b*c)))/2
# set up result screen
axis((0,10,0,10))
axis("off")
# determine solutions
# best used for real solutions
if p!=q:
r=(b*n+m*(p-d))/(p-q)
s=m-r
u=(c*m+n*(p-a))/(p-q)
v=n-u
text(0,9,"x=re**(pt)+se**(qt)")
text(0,8.5,"r="+str(r))
text(0,8,"p="+str(p))
text(0,7.5,"s="+str(s))
text(0,7,"q="+str(q))
text(0,6,"y=ue**(pt)+se**(vt)")
text(0,5.5,"u="+str(u))
text(0,5,"p="+str(p))
text(0,4.5,"v="+str(v))
text(0,4,"q="+str(q))
else:
r=(b*n+m*(p-d))
s=m
u=(c*m+n*(p-a))
v=n
text(0,9,"x=rte**(pt)+se**(pt)")
text(0,8.5,"r="+str(r))
text(0,8,"p="+str(p))
text(0,7.5,"s="+str(s))
text(0,6,"y=ute**(pt)+se**(vt)")
text(0,5.5,"u="+str(u))
text(0,5,"p="+str(p))
text(0,4.5,"v="+str(v))
show()
Examples
Example 1
dx/dt = 3 * x(t) - 2 * y(t)
dy/dt = 2 * x(t) - 2 * y(t)
x(0) = 6, y(0) = 9
Solutions:
x(t) = 2 * e^(2*t) + 4 * e^(-t)
y(t) = e^(2*t) + 8 * e^(-t)
Example 2
dx/dt = 3 * x(t) - 4 * y(t)
dy/dt = x(t) - y(t)
x(0) = 1, y(0) = 0
Solutions:
x(t) = 2 * t * e^(t) + e^(t)
y(t) = t * e^(t)
Eddie
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.