TI-84 Plus CE and Numworks: System of Linear Differential Equations (Simple System)

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.