===================================================
Math 451
Laboratory 4: Differential Equations
Date: 2013-04-05
Due: 2013-04-12
===================================================
Differential equations are one of the most useful conceptual
tools of applied mathematics. But solving them requires
some practice.
Goals
----------------------------
- Learn how to solve a differential equation with Matlab
- Learn to nth order equations as systems
- Learn to plot time series and phase portraits
- Solve Newton's equations for a modified 3-body problem
- Make javascript animations of differential equation
- Solve Newton's planet equations with python
Note: The differential equation solver in Octave (lsode) is
not compatible with Matlab, so you will have to do some
creative re-writing as you go if you would like to complete
this lab using Octave.
===================================================
A simple example of an initial value problem
----------------------------------------------
The first-order equation
dn/dt = r n (1 -n ) - k (1 - cos(2*pi*t) )
is what's called a Ricatti equation. The variable n(t) represents the
density of oysters relative to their natural carrying capacity and the
parameter r represents the maximum per-capital growth rate of the
population. The oysters compete for space on the ocean bottom, and
they are harvested at different rates over the course of the year. We
would like to know if harvesting rates corresponding to k = 0.2 are
too intense and will collapse the stock. To do this, we'd like to
explore the solutions of this equation, but unfortunately, this
equation does not have closed-form solutions. So instead, we will
solve it numerically.
To solve this equation numerically, we need to specify the equation's
right-hand side, an initial condition, and a sequence of time-points
where we would like to calculation the value of the solution curve.
(1-a) Create a Matlab function file called "oyster_ode.m" This
function takes the current time and state as arguments and returns
a vector representing the direction and rate of change in the next
time-step.
function dn = oyster_ode(t,n)
k = 0.2;
r = 1.0;
dn = [r*n(1)*(1 - n(1))-k*(1-cos(2*pi*t))];
return
(1-b) Create a Matlab script called 'oyster_simulator.m' using the
code below.
ic = [ 1. ];
tspan = linspace(0,8,1000);
[t,x] = ode45( 'oyster_ode', tspan, ic);
figure(1);
plot(t,x,'o-');
xlabel('Time');
ylabel('Oyster density');
ylim([-1,3]);
This file specifies and solves the differential equation. Here, ic
sets the initial oyster population size and tspan is the initial and
final times at which we would like to know solution. ode45 is a
numerical method for solving differential equations. Specifically, it
is a 5th order Runge-Kutta method with adaptive time stepping. It
returns a vector t of times, and a matrix x of the solution at those
times.
(1-c) Run your script 'oyster_simulator.m'.
(1-d) Does the stock crash and go extinct?
(1-e) To show how what happens depends on the initial population size,
plot solutions for 20 initial conditions ranging from 0 to 3 on
the same plot. Save your figure to be handed in.
You may find the "hold on" command useful.
Solutions of systems
--------------------------------------------------------
Numerical solvers like ode45 can only solve 1st-order equations. So
whenever we need to solve a 2nd order equation, we transform it into a
system of first order equations and solve those.
Consider the damped pendulum equation
u" = - c u' - sin( u )
with c = 0.1. This is a second order nonlinear equation.
(2-a) Use the substitution y(1) = u, y(2) = u' to rewrite this
equation as a system of 2 first-order equations.
(2-b) Write a new function file called "pendulum_vfield.m"
in the same style as "oyster_ode.m" that takes a time and a
position y as input arguments and returns a vector representing
the current state velocity AS A COLUMN VECTOR.
(2-c) Write a new script named "pendulum_simulator.m" that
plots time-series for u(t) and u'(t) starting from an initial
position u(0) = pi/2 and initial velocity u'(0) = 0 over 25
seconds. Include labels and a legend. Save this figure to hand in.
Phase portraits are plots of the states of 2-dimensional systems
relative to each other. For some simple examples, see
http://tutorial.math.lamar.edu/Classes/DE/PhasePlane.aspx
and http://en.wikipedia.org/wiki/Phase_portrait
(2-d) Construct solutions to the pendulum equation with u'(0)
is in linspace(-5,5,15) and u(0) = 0, and plot these solutions
to make a phase portrait. Label your axes and save this
figure to hand in.
(2-e) Why are the spirals repeating?
(optional) What are the equilibrium points and where are they in the
phase portrait? What are the basins of attraction of attracting
equilibria?
Javascript and Differential Equations in Browsers
----------------------------------------------------
One downside of the images we've produced so far is that they are
static, and don't communicate any of the dynamic nature of
differential equations. But our web browsers with javascript and
HTML5 are now very powerful tools for doing this.
For example, see http://hint.fm/wind/
Let's do a simpler version of this (the original version of this is
at https://gist.github.com/1320178).
(3-a) Download the "lorenz_equation.html" to your Desktop, open it,
and read it. The file contains a combination of html markup,
style-sheet specifications, and javascript code.
* The first part in the "style" tags is the style sheet stuff
* The next part is the script to solve our differential equation.
Even though we haven't used javascript in this class before,
you should be able to recognize variable names and equations
to get an idea for how the code works.
(3-b) Reconstruct the specification of the differential equations
this code is solving using Euler's method.
(3-c) What is the step size being used in Euler's method?
Now, test the code in a web browser to see what it does.
(3-d) You may notice that the dots change size. What does this mean?
(3-e) How many particles are being used in the paths?
Now, we'll make some simple changes to the code to change the
animations we get out.
(3-f) How does the behavior of the animation change as you use large
step sizes? Be specific. (Don't increase the step size too fast,
or it will be hard to tell.)
(3-g) Modify the code so that we can see the cumulative shape of the
attractor. Describe your changes, and include a screen-shot of
the result.
Newton's equations for planetary motion in python
----------------------------------------------------
The original killer application of differential equations
was explaining planetary motion. Unfortunately, these
equations are non-linear and at least 4-dimensional. This
is too hard for us to solve in our first differential
equations class, but they are not too hard to solve
numerically! The main concept is that forces from gravity
on a body are proportional to the masses of both bodies and
inversely proportional to the distances between the bodies.
If r is the radial position vector of a planet from it's
sun, then Newton's equation is
2
d r G m r
---- = -------
2 3
dt |r|
(4-a) Rewrite this equation as a system of 4 1st-order
differential equations where r = (x,y).
We will solve these equations using a python script.
(4-b) Download "planetary_motion.py".
This script predicts the motional of a small satellite
in a system with a star (or starts) that have fixed
positions.
(4-c) Run the script "planetary_motion.py" in the terminal
using the command
python planetary_motion.py
One figure you get should be the path of a planet around
1 sun. The other is the time-series for the x and y
positions and the x and y velocities. Does the picture
agree with your and Johannes Kepler's expectations?
(4-d) Add a second sun to the system, with the same mass as
the first sun, but at position [-0.5,0.]. Keep everything
else the same.
(4-e) Run your new script and save the orbital path figure to
hand in.
(Extra credit) Port this script to Matlab.