Mathematica  for Systems of Linear Differential Equations

Mathematica can solve systems of differential equations. As always, some analysis of the solution which Mathematica returns may be in order.

Homogeneous System

deq1 = x1'[t] ==-2x1[t] + x2[t]
deq2 = x2'[t]==x1[t]-2x2[t]

DSolve[{deq1,deq2},{x1[t],x2[t]},t]

Nonhomogeneous System

deq1 = x1'[t] ==-2x1[t] + x2[t] + 2 Exp[-t]
deq2 = x2'[t]==x1[t]-2x2[t] + 3 t

DSolve[{deq1,deq2},{x1[t],x2[t]},t]

Initial Value Problem System

deq1 = x1'[t] ==-2x1[t] + x2[t]
deq2 = x2'[t]==x1[t]-2x2[t]

DSolve[{deq1,deq2,x1[0]==3/4,x2[0]==-1},{x1[t],x2[t]},t]

Eigenvalues and Eigenvectors

As always, it is better to know something behind what you are doing than just having Mathematica do it for you (although you may want to use Mathematica to check your results).

We formulated solutions to systems of linear differential equations based on the eigenvalues and eigenfunctions of the transformation matrix A. For example, the above homogeneous problem could be written as:

x[t_]= {x1[t],x2[t]}
A = {{-2,1},{1,-2}}
MatrixForm[x[t]]
MatrixForm[A]
deq = D[x[t],t] == Dot[A ,x[t]]

We assume a solution looks like

Solution[t] = {ζ1,ζ2} Exp[λ t]

where the vector {ζ1,ζ2} is the eigenvector associated with the eigenvector λ for the matrix A. We can have Mathematica find these for us:

Eigenvectors[A]
Eigenvalues[A]
Eigensystem[A]

So the general solution looks like:

sol1 = c1 {-1,1} Exp[-3 t] + c2{1,1} Exp[-1 t]

This is the same as found above, although at first glance they may not look the same.

Creating Direction Fields for Systems

A = {{-2,2},{1,-2}}

field=VectorPlot[A.{x1,x2},{x1,-3,3},{x2,-3,3},ImageSize->{500,Automatic},FrameLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},VectorScale->{Tiny,Small,None},VectorPoints->30,BaseStyle->20]

FYI: Mathematica can deal with the exponential of a matrix!

Let's revisit an earlier example, and see what happens if we try to work with the exponential of a matrix. This is strictly for fun, since this is not how we will be solving these problems in this class.

Here is the Mathematica solution to a linear system initial value problem:

deq1 = x1'[t] ==5x1[t] -x2[t]
deq2 = x2'[t]==3x1[t]+x2[t]
DSolve[{deq1,deq2,x1[0]==-1,x2[0]==2/3},{x1[t],x2[t]},t]

Here is another solution, based on the exponential of a matrix. We identify the matrix A:

A = {{5,-1},{3,1}}

The solution is given by:

sol[t_]=MatrixExp[A t].{c1,c2}
Collect[%,{c1,c2}]

Use the initial conditions to determine the constants:

constants=Solve[sol[0]=={-1,2/3},{c1,c2}]

The solution to the IVP is:

sol[t_]=Simplify[sol[t]/.Flatten[constants]]

which is the same as Mathematica returned! So it might be a useful thing in the future to understand the exponential of a matriax, but that requires a bit more linear algebra than this course requires.

ParametricPlot[sol[t],{t,-1,1},Frame-> True,ImageSize->{500,Automatic},FrameLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},BaseStyle->20]