C Documents and Settings Administrator Ustawienia lokalne Dane aplikacji Opera Opera cache g 0021 opr00FZT

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 y

′

(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 y

′

= y. A solution is given by y(t) = e

(just plug it in and check to

see).

More generally, any function of the form y(t) = ce

satisfies y

′

= y; we say that y = ce

is the

general solution

to y

′

= y. We can specify the solution precisely using an initial condition. For

example, y(0) = 2 together with y

′

= y implies that y(t) = 2e

. To find this, we would first write

(t) = ce

and then use the initial condition y(0) = 2 to find that c = 2. The pair of equations

′

= 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 commands Solve (a purely algebraic, exact
solver) and NSolve (an approximate, numeric solver). There are analogs of these for differential
equations, namely DSolve and NDSolve.

DSolve

attempts to find an exact solution of an ordinary differential equation or a system of

ordinary differential equations. The syntax is DSolve[eqs, funcs, var], where eqs is the list of
equations to solve for, funcs is the unknown function or functions, and var is the variable of the
unknown function. Here are few examples illustrating this syntax.

The most basic example:

DSolve

y '

D 

, y

, t

We can add an initial condition.

DSolve

y '

D 

, y

D 

, y

, 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 function y[t] as well. It’s usually a good idea to Clear the function
first to avoid any problems that might arise, if it’s already defined.

Clear

;

sol

< =

DSolve

y '

D 

D H

, y

D 

, y

, t

;

D =

D 

. 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 in y[t]. Just plug in to check.

Expand

y '

, y

D H

DL<D

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 y

′

= f (y, t). Suppose that the graph of the solution y(t) passes through

the point (t

, y

). Then it must do so with slope y

′

), which must be f (y

, t

) 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 the VectorFieldPlot command from the
VectorFieldPlots

package.

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 the Needs command that loads the package and
replace the VectorFieldPlot command with either VectorPlot or StreamPlot. The result will
be much nicer, particularly if you use StreamPlot.

At any rate, here’s the V6 and V7 compatible code to generate a plot of a vector field.

Needs

"VectorFieldPlots`"

;

VectorFieldPlot

1, y

4, 4

y, 0, 4

Now, let’s show the vector field together with a graph of the solution.

Plot

3, 3

PlotStyle

Thick

;

Show

g, vf

Note how the solution simply follows the arrows.

Exercise 1:

Use Mathematica to solve the IVP y

′

= 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,

+ y

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.

′

(t)

= 2x(t) − y(t)

′

(t)

= x(t) + 2y(t)

We can solve systems using DSolve as 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 '

D 

D -

y '

D 

D +

D 

1, y

D -

, y

, 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 (x

′

(t), y

′

(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.

, y

D< = 8

, y

D< 

. sol;

ParametricPlot

, y

2, 2

PlotStyle

Thick

;

VectorFieldPlot

y, x

4, 4

;

Show

g, vf

, PlotRange

88-

4, 4

<<D

This example looks really cool with V7’s StreamPlot.

Exercise 2:

Solve the system

′

(t)

= 2x(t) − y(t)

′

(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 solvers NSolve and FindRoot for 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 solver NDSolve is

very broadly applicable, however.

We start by looking at our most basic example again: y

′

= y, y(0) = 2. We know that this IVP

has the solution y(t) = 2e

, so it will be easy to check the accuracy of our numerical approximator.

Here’s the syntax for solving this IVP numerically.

Clear

;

sol

< =

NDSolve

y '

D 

, y

D 

, y

2, 2

A few comments are in order here. First, the syntax is very similar to the DSolve syntax for solving
an IVP. The only difference is that the independent variable specification changes from a simple t
for DSolve to {t,-2,2} for NDSolve. This specifies the domain on which our approximation will
be valid. Second, the result looks very strange indeed. It’s an InterpolatingFunction object,
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, an InterpolatingFunction can 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 function y[t].

D =

D 

. sol

Now, we can compute the value at any point in [−2, 2]. The value at 0 should be 2, right?

Cool. Here’s the value at 2.

This should be 2e

Looks good so far. Let’s try to plot it.

Plot

2, 2

Looks just like the graph of an exponential function! The point is that we can treat y just 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.

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 2e

, you’ll

find that it’s quite a bit off.

Exercise 3:

Use NDSolve to plot the solution to the IVP y

′

= y sin y

, 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

′

(y − x)

′

= ρx − xz − y

′

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 using NDSolve as follows, though.

Clear

x, y, z

;

sol

< =

NDSolve

x '

D 

D -

y '

D 

D -

z '

D 

D -

D 

0, y

D 

1, z

D 

, y

, z

t, 0, 100

;

, y

, z

D< = 8

, y

, z

D< 

. sol;

ParametricPlot3D

, y

, z

t, 0, 100

Now, that’s cool!

Exercise 4:

The R¨

ossler attractor is the orbit of the solution of the non-linear system

′

−(y + z)

′

+ ay

′

+ z(x − c).

Plot the solution of this system using the parameters a = 0.2, b = 0.2 and c = 5.7.

Wyszukiwarka

Podobne podstrony:

więcej podobnych podstron