Intro to Mathematica - Lab 10
Differential equations
Due to the wide variety of backgrounds in this class, we ve largely avoided specific topics of study.
The reality is, however, that Mathematica is primarily a tool for doing math and differential
equations is a topic that virtually everyone will see that benefits greatly from computational power.
While the DE class is not a prerequisite for this class, the hard work of solving the equations can
be done by Mathematica and the basic concepts are easily explained.
Note that differential equations is suffused with acronyms - DE for differential equation and IVP
for initial value problem, for example. You ll get the hang of it.
The basics
Some definitions
A first order, ordinary differential equation has the form y2 (t) = f(y, t). A solution to this differ-
ential equation is a function y(t) that satisfies the equation. The most basic example is given by
f(y, t) = y; thus the DE is y2 = y. A solution is given by y(t) = et (just plug it in and check to
see).
More generally, any function of the form y(t) = cet satisfies y2 = y; we say that y = cet is the
general solution to y2 = y. We can specify the solution precisely using an initial condition. For
example, y(0) = 2 together with y2 = y implies that y(t) = 2et. To find this, we would first write
y(t) = cet and then use the initial condition y(0) = 2 to find that c = 2. The pair of equations
y2 = y, y(0) = 2 is called an initial value problem or IVP.
While differential equations can be much richer and more complicated than this simple example
suggests, these basic definitions are enough to get started.
DSolve
When we solve regular equations, we frequently use the commandsSolve(a purely algebraic, exact
solver) andNSolve(an approximate, numeric solver). There are analogs of these for differential
equations, namelyDSolveandNDSolve.
DSolveattempts to find an exact solution of an ordinary differential equation or a system of
ordinary differential equations. The syntax isDSolve[eqs,funcs, var], whereeqsis the list of
equations to solve for,funcsis the unknown function or functions, andvaris the variable of the
unknown function. Here are few examples illustrating this syntax.
The most basic example:
DSolve y' t y t , y t , t
We can add an initial condition.
DSolve y' t y t , y 0 2 , y t , t
Here s a more complicated example. In this example, we solve for the function y(t) and we store
the result in the Mathematica functiony[t]as well. It s usually a good idea toClearthe function
first to avoid any problems that might arise, if it s already defined.
Clear y ;
sol DSolve y' t y t 4 y t , y 0 2 , y t , t ;
y t y t . sol
Note that a warning message is generated, so it s probably not a bad idea to check the solution.
This is easy to do, since we ve got the solution iny[t]. Just plug in to check.
Expand y' t , y t 4 y t
Looks like they re the same.
Visualizing the solution
We can visualize the solution to a differential equation by simply graphing the solution as we
would any function. The solution is more understandable, if we graph it on top of the slope field
generated by the differential equation. To understand this, consider the geometric interpretation
of the differential equation y2 = f(y, t). Suppose that the graph of the solution y(t) passes through
the point (t0, y0). Then it must do so with slope y2 (t0), which must be f (y0, t0) by the differential
equation. To visualize this, we simply place an arrow with slope f(y, t) at each point (t, y) chosen
from a grid of points. This can be done easily using theVectorFieldPlotcommand from the
VectorFieldPlotspackage.
Note: Mathematica V7 introduces greatly improved functions for plotting vector fields. The fol-
lowing code is V6 compatible and works in V7 as well, but V7 will warn you that there is a better
way. If you are using V7, you can simply delete theNeedscommand that loads the package and
replace theVectorFieldPlotcommand with eitherVectorPlotorStreamPlot. The result will
be much nicer, particularly if you useStreamPlot.
At any rate, here s the V6 and V7 compatible code to generate a plot of a vector field.
Needs "VectorFieldPlots`" ;
vf VectorFieldPlot 1, y 4 y , t, 4, 4 , y, 0, 4
Now, let s show the vector field together with a graph of the solution.
g Plot y t , t, 3, 3 ,
PlotStyle Thick ;
Show g, vf
Note how the solution simply follows the arrows.
Exercise 1: Use Mathematica to solve the IVP y2 = e-t - y and plot the solution on top of the
slope field.
Systems of differential equations
In elementary algebra, you learn about systems of equations. For example,
x + y = 2
x - y = 0
is a linear system of two equations in the two unknowns x and y. The unique solution is x = 1,
y = 1.
In a similar manner, we can work with systems of differential equations. Here s an example of a
first order, linear system of differential equations.
x2 (t) = 2x(t) - y(t)
y2 (t) = x(t) + 2y(t)
We can solve systems usingDSolveas easily as we can solve single equations. Here s the syntax
to solve this system subject to the initial condition x(0) = 1, y(0) = -1.
Clear x, y ;
sol DSolve
x' t x t y t ,
y' t x t y t ,
x 0 1, y 0 1 ,
x t , y t , t
As before, we d like to visualize this solution on top of a relevant vector field. Now, the solution
may be thought of as a parametric function p(t) = (x(t), y(t)); as t moves, the point p(t) traces out
a curve in the plane. The velocity of the point is (x2 (t), y2 (t)) which, by the differential equations
is (x - y, x + y). Thus the motion must follow the vector field F (x, y) = (x - y, x + y). We can
illustrate this as follows.
x t , y t x t , y t . sol;
g ParametricPlot x t , y t , t, 2, 2 ,
PlotStyle Thick ;
vf VectorFieldPlot x y, x y , x, 4, 4 , y, 4, 4 ;
Show g, vf , PlotRange 4, 4 , 4, 4
This example looks really cool with V7 sStreamPlot.
Exercise 2: Solve the system
x2 (t) = 2x(t) - y(t)
y2 (t) = x(t) + 2y(t)
subject to the initial condition x(0) = 1, y(0) = -1 and plot the solution on top of the relevant
vector field.
NDSolve
Recall that many simply stated equations can t be solved exactly. Mathematica provides the
numerical solversNSolveandFindRootfor such situations. Similarly, not all differential equations
can be solved exactly. In fact, important differential equations arising from applied situations
usually cannot be solved in closed form. The numerical differential equation solverNDSolveis
very broadly applicable, however.
We start by looking at our most basic example again: y2 = y, y(0) = 2. We know that this IVP
has the solution y(t) = 2et, so it will be easy to check the accuracy of our numerical approximator.
Here s the syntax for solving this IVP numerically.
Clear y ;
sol NDSolve y' t y t , y 0 2 , y t , t, 2, 2
A few comments are in order here. First, the syntax is very similar to theDSolvesyntax for solving
an IVP. The only difference is that the independent variable specification changes from a simplet
forDSolveto{t,-2,2}forNDSolve. This specifies the domain on which our approximation will
be valid. Second, the result looks very strange indeed. It s anInterpolatingFunctionobject,
which simply stands for a numerical procedure to approximate the solution at any point in the
domain. You learn about one such procedure (called Euler s method) in a differential equations
class, but the details are not particular important for us. On the other hand, we do need to know
how to access the solution. As turns out, anInterpolatingFunctioncan be used just as any
other function, even though we have no simple formula for it. For example, let s store the result
in the functiony[t].
y t y t . sol
Now, we can compute the value at any point in [-2, 2]. The value at 0 should be 2, right?
y 0
Cool. Here s the value at 2.
y 2
This should be 2e2.
N 2 2
Looks good so far. Let s try to plot it.
Plot y t , t, 2, 2
Looks just like the graph of an exponential function! The point is that we can treatyjust as we
would any function in Mathematica; we just don t have a formula for it. We do need to be careful
to stay inside our specified domain of definition, though. Here s what happens otherwise.
y 5
Note that a value is computed by extrapolation, but a warning is issued and there s no reason to
expect the computed value to be very good. If you check against the actual value of 2e5, you ll
find that it s quite a bit off.

Exercise 3: UseNDSolveto plot the solution to the IVP y2 = y sin y2 , y(0) = 1.
Non-linear systems
Non-linear systems of differential equations have been an active field of research for literally cen-
turies. They arise already in the work of Isaac Newton on celestial mechanics, yet their study
continues to be a growing field. Perhaps the most famous non-linear system from the last half-
century is the Lorenz system. This arises in the field of atmospheric dynamics and has become an
icon of chaos. The Lorenz system is
x2 = Ã(y - x)
y2 = Áx - xz - y
z2 = xy - ²z.
The symbols Ã, Á, and ² are constant parameters. The symbols x, y, and z are functions of the
independent variable t. Thus, the solution p(t) = (x(t), y(t), z(t)) traces out a path in three-space.
There is no closed form for the solution and, for certain values of the parameters, the solution is
known to be chaotic. We can approximate the solution usingNDSolveas follows, though.
Clear x, y, z ;
sol NDSolve
x' t 3 y t x t ,
y' t 26 x t x t z t y t ,
z' t x t y t z t ,
x 0 0, y 0 1, z 0 1 ,
x t , y t , z t ,
t, 0, 100 ;
x t , y t , z t x t , y t , z t . sol;
ParametricPlot3D x t , y t , z t , t, 0, 100
Now, that s cool!
Exercise 4: The Rössler attractor is the orbit of the solution of the non-linear system
x2 = -(y + z)
y2 = x + ay
z2 = b + z(x - c).
Plot the solution of this system using the parameters a = 0.2, b = 0.2 and c = 5.7.