MATH 302-010
Ordinary Differential Equations
Prof. D. A. Edwards
Copyright 2024

Advanced Linear Systems

We begin with our standard command:

Quit[]

Since we are working with matrices and vectors, we also insert the following command, which makes them display as they do in class:

$PrePrint=If[MatrixQ[#] || VectorQ[#],MatrixForm[#],#]&;

Complex Eigenvalues

All the results for real eigenvalues will work for complex eigenvalues, where we recall that Mathematica denotes (Sqrt[-1]) by ι:

b={{2,1},{-1,3}}
Eigensystem[b]
bval = %[[1]]
bvecs=%%[[2]]

To check that we do indeed get the right eigenvectors, we calculate Bz to see if it is indeed  λ z:

lam1 = bval[[1]]
z1 = bvecs[[1]]
ComplexExpand[b.z1]
ComplexExpand[lam1*z1]

Exercise

1.  Find the determinant, eigenvalues, and eigenvectors of the following matrix:

A=(({{0, 4}, {-13, 9}}))

Solution

mat1 = {{0,4},{-13,9}}
Det[mat1]
Eigensystem[mat1]

Phase Plane (Complex Eigenvalues)

Solving systems of ODEs with oscillatory solutions works much as before.

Let’s consider the following system with complex eigenvalues:

mat2= {{3,-2},{4,2}}
Eigensystem[mat2]

Again we should use ComplexExpand to properly interpret our expressions:

ComplexExpand[%]

If you have a matrix with complex eigenvalues, Mathematica still knows that the solutions of the ODE should be real:

sys2a=x'[t]==3*x[t]-2*y[t]
sys2b=y'[t]==4*x[t]+2*y[t]
DSolve[{sys2a,sys2b},{x[t],y[t]},t]
ComplexExpand[%]

Similarly, the StreamPlot command produces a real direction field:

sp2={3*x-2*y,4*x+2*y}
StreamPlot[sp2,{x,-1,1},{y,-1,1}]

Note that we have an unstable spiral, as suggested by the eigenvalues.

Exercises

1.  Find the general solution of the following system:

x'=x+5y
y'=-2x+7y

Then graph the phase plane with StreamPlot.

Solution

We start by creating the form that StreamPlot needs:

p1r1 = x+5*y
p1r2 = -2*x+7*y
StreamPlot[{p1r1,p1r2},{x,-2,2},{y,-2,2}]

Then we put our system into the form that DSolve needs.

p1eq1 = x'[t]==x[t]+5*y[t]
p1eq2 = y'[t]==-2*x[t]+7*y[t]

DSolve[{p1eq1,p1eq2},{x[t],y[t]},t]

Repeated Eigenvalues

Solving systems of ODEs with repeated solutions works much the same as before.

If a matrix has repeated eigenvalues, Mathematica lists (0,0) as the “second” eigenvector:

matr= {{1,-2},{2,-3}}
Eigensystem[matr]

However, Mathematica still knows the proper form of the solutions:

xvec={x[t],y[t]}
sys1a = x'[t]==x[t]-2*y[t]
sys1b=y'[t]==2*x[t]-3*y[t]
DSolve[{sys1a,sys1b},xvec,t]

Similarly, the StreamPlot command produces a real direction field:

sp1={x-2*y,2*x-3*y}
StreamPlot[sp1,{x,-1,1},{y,-1,1}]

Note that we have only a single eigenvector.

Exercise

1.  Find the general solution of the following system:

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

Then graph the phase plane with StreamPlot.

Solution

We write the system in the form that StreamPlot needs:

p1r1 = x+2*y
p1r2 = -2*x-3*y
StreamPlot[{p1r1,p1r2},{x,-2,2},{y,-2,2}]

Then we put our system into the form that DSolve needs.

p1eq1 = x'[t]==x[t]+2*y[t]
p1eq2 = y'[t]==-2*x[t]-3*y[t]

DSolve[{p1eq1,p1eq2},{x[t],y[t]},t]

Inhomogeneous Systems of Linear ODEs

Mathematica can also solves systems of forced ODEs.  In particular, consider the following set of equations:

fsys1a=x'[t]==3*x[t]+2*y[t]+Exp[-2*t]
fsys1b = y'[t]==-2*x[t]-2*y[t]-Sin[t]

Again we use DSolve, listing both variables to be found as a set, as shown below:

DSolve[{fsys1a,fsys1b},{x[t],y[t]},t];
Simplify[%]

Note that two arbitrary constants appear, but they appear in BOTH functions, since they are related to one another by the system.  If you have particular initial conditions to solve, you can add them---as many or as few as you have:

ic1a=x[0]==3
ic1b=y[0]==5
DSolve[{fsys1a,fsys1b,ic1a,ic1b},{x[t],y[t]},t]

Real or complex eigenvalues work the same:

mat2= {{3,-2},{4,2}}
Eigensystem[mat2]

fsys2a=x'[t]==3*x[t]-2*y[t]+Exp[2*t]
fsys2b = y'[t]==4*x[t]+2*y[t]-Sin[t]
DSolve[{fsys2a,fsys2b,ic1a,ic1b},{x[t],y[t]},t];
Simplify[%]

Exercise

1.  Find the general solution of the following system:

x'=x-y+(e^(-t))
y'=-2x+7y+sin(3t)

Then find the solution subject to x(0)=3, y(0)=5.

Solution

p1eq1 = x'[t]==x[t]-y[t]+Exp[-t]
p1eq2 = y'[t]== -2*x[t]+7*y[t]+Sin[3*t]
DSolve[{p1eq1,p1eq2},{x[t],y[t]},t]

p1ic1 = x[0]==3
p1ic2 = y[0]==5
DSolve[{p1eq1,p1eq2,p1ic1,p1ic2},{x[t],y[t]},t]

Nonlinear Systems of ODEs

Mathematica is not as good at solving nonlinear sets of equations.  (When you execute this line, Mathematica will seem to hang—just hit CMD-. to stop it.)

nlsys1 = x'[t]==-4*(x[t])^2-3*y[t]
nlsys2 = y'[t]==3*Sin[x[t]]-2*(y[t])^3
DSolve[{nlsys1,nlsys2},xvec,t]

Thus we use the NDSolve option, though of course we have to give initial conditions:

nlsyssol = NDSolve[{nlsys1,nlsys2,x[0]==1,y[0]==-1},{x,y},{t,0,1}];

Again note that we drop the arguments of x and y.

Then as before, if we put in a particular value of t, we can get the solution for BOTH x and y:

xvec/.(t->1)
%/.nlsyssol

We can also create phase-plane plots.  Here we widen the view somewhat to capture both fixed points.

spnl1={-4*x^2-3*y,3*Sin[x]-2y^3}
StreamPlot[spnl1,{x,-2,2},{y,-2,2}]

Mathematica can solve forced nonlinear sets of equations too, though now we are wise enough to proceed immediately to the numerical form:

fnlsys1=x'[t]==-4*(x[t])^2-3*y[t]+Sin[t]
fnlsys2=y'[t]==3*Sin[x[t]]-2*(y[t])^3-2*Exp[-t]
fnlsyssol= NDSolve[{fnlsys1,fnlsys2,x[0]==1,y[0]==-1},{x,y},{t,0,1}];

In the forced case, phase planes are not particularly helpful.  But we can plot one function at a time:

Plot[x[t]/.fnlsyssol,{t,0,1},AxesLabel->{t,x}]
Plot[y[t]/.fnlsyssol,{t,0,1},AxesLabel->{t,y}]

Exercises

1.  On the same graph, plot the solutions x(t) and y(t) vs. t for the following system for 0<t<5:

x'=sin(x) - 3y
y'=cos(x)-2(y^3)
x(0)=0, y(0)=0

Solution

p2eq1 = x'[t]==Sin[x[t]]-3*y[t]
p2eq2 = y'[t]==Cos[x[t]]-2*(y[t])^3
p2ic1 = x[0]==0
p2ic2 = y[0]==0
p2sol = NDSolve[{p2eq1,p2eq2,p2ic1,p2ic2},{x,y},{t,0,5}];

plex1a=Plot[x[t]/.p2sol,{t,0,5}]
plex1b=Plot[y[t]/.p2sol,{t,0,5}]
Show[{plex1a,plex1b},PlotRange->All]

2.  Graph the phase plane of the following system:

x'=sin(x) - 3y
y'=cos(x)-2(y^3)

Solution

We write the system in the form that StreamPlot needs:

p2r1 = Sin[x]-3*y
p2r2 = Cos[x]-2*y^3
StreamPlot[{p2r1,p2r2},{x,-2,2},{y,-2,2}]

3.  On the same graph, plot the solutions x(t) and y(t) vs. t for the following system for 0<t<5:

x'=sin(x) - 3y+sin((t^3))
y'=cos(x)-2(y^3)-(e^(-(t^2)))
x(0)=0, y(0)=0

Solution

p2eq1 = x'[t]==Sin[x[t]]-3*y[t]+Sin[t^3]
p2eq2 = y'[t]==Cos[x[t]]-2*(y[t])^3-Exp[-t^2]
p2ic1 = x[0]==0
p2ic2 = y[0]==0
p2sol = NDSolve[{p2eq1,p2eq2,p2ic1,p2ic2},{x,y},{t,0,5}];

plex3a=Plot[x[t]/.p2sol,{t,0,5}]
plex3b=Plot[y[t]/.p2sol,{t,0,5}]
Show[plex3a,plex3b,PlotRange->All]

Laplace Transforms

To use Mathematica for Laplace transforms, we use the command LaplaceTransform[f(t),t,s], where t is the regular variable in regular space and s is the transformed variable:

f1 = f[t]==Exp[-t]
eq1 = LaplaceTransform[f1,t,s]

Note that Mathematica indicates the transform (which I call (OverHat[f])) by (Subscript[ℒ, t]) (depending on how your formatting is set up).    It is important to put the variables in the proper order, or you will get back garbage:

LaplaceTransform[f1,s,t]

Mathematica can’t do all possible Laplace transforms, and will simply return the input if the transform is too complicated:

f2 = Exp[-Sin[t]]
LaplaceTransform[f2,t,s]

To invert the Laplace transform, we use the command InverseLaplaceTransform[fhat(s),s,t]:

InverseLaplaceTransform[eq1,s,t]

Note that Mathematica can invert even symbolic transforms like the left-hand side of the above.  Now consider a complicated fractional transform:

eq2=1/s/(s+1)/(s^2+4)/(s^4-1)
InverseLaplaceTransform[eq2,s,t]

To see where the parts come from (or to expand such an expression to invert by hand), use the Apart command, which does partial fractions:

Apart[eq2,s]

Exercises

1.  Calculate the Laplace transform of (e^(-3t))+(e^(-t))sin(t).

Solution

prob1 = Exp[-3*t]+Exp[-t]*Sin[t]
LaplaceTransform[prob1,t,s]

2.  Calculate the inverse Laplace transform of (1/(s+3))+(3/((s^2)+4s)).

Solution

prob2 = 1/(s+3)+3/(s^2+4*s)
InverseLaplaceTransform[prob2,s,t]