Mathematical Modeling: laboratory exercises

Online documentation for software

Online help for matlab is available in, e.g.,


gnuplot is a simple to use (and free) software to plot graphs, contours, etc. Online help is available in, e.g.,

Modeling software and online examples

Population models

Exercise (Lotka-Volterra model)
Study the Lotka-Volterra model.
(a) Find the equilibrium points for a,b,c,d of your choice (e.g., a=3, b=10, c=1, d=2).
(b) Solve numerically the system of equations for various initial conditions (use, e.g., matlab).
(c) Plot the populations as a function of time.
(d) Plot the phase portrait.

(a) See the Lecture Notes.

(b), (c) Here is a guide (from MathWorks website).
You can find another guide in the book [Barnes & Fulford], page 109.

Here is a basic code in matlab:

global a b c d;
a=10; b=2; c=3; d=1;
t_init=0;
t_fin=5;
hold on;
y0 = [2 3];
[t,y] = ode45('lot_volt',[t_init,t_fin],y0);
plot(t,y(:,1));
hold on
plot(t,y(:,2));
axis([0 5 0 10]);
xlabel('time');
ylabel('x, y');
title('Lotka-Voltera solution');

In the above code we use the function:

function yp=lot_volt(t,y)
global a b c d;
h1 = a*y(1) - c*y(1)*y(2);
h2 = -b*y(2) + d*y(1)*y(2);
yp = [h1 h2]';

Here is the result:
Lotka-Voltera solution

(d) We should use the integral for the model given in the lecture notes and make a contour plot for it. We can use gnuplot and similar commands is in Exercise "Simple phase portrait".

Alternatively, we may use matlab to directly solve the differential equations (and then plot the phase portrait), for example:

global a b c d;
a=10; b=2; c=3; d=1;

t_init=0;
t_fin=1;
hold on;
y0 = [5 5]; [t,y] = ode45('lot_volt',[t_init,t_fin],y0);
plot(y(:,1),y(:,2));
y0 = [6 5]; [t,y] = ode45('lot_volt',[t_init,t_fin],y0);
plot(y(:,1),y(:,2));
y0 = [7 5]; [t,y] = ode45('lot_volt',[t_init,t_fin],y0);
plot(y(:,1),y(:,2));
y0 = [8 5]; [t,y] = ode45('lot_volt',[t_init,t_fin],y0);
plot(y(:,1),y(:,2));
axis([0 9 0 9]);
xlabel('x');
ylabel('y');
title('Lotka-Voltera phase diagram');

Here is the result:
Lotka-Voltera phase portrait

Phase portraits

Exercise (Simple phase portrait)
Plot the phase portrait for
dx/dt = ay, dy/dt= -x

Fixed Point: (i) x=0, y=0

Phase curves are given by
x2/2+a y2/2 = C (for all C>=0)
[explain why.]

Let us first plot the function (choose. e.g., a=2) f(x,y)= x2/2+a y2/2.
Using gnuplot we have the commands:

a=2.0
set xrange[-10:10]
set yrange[-10:10]
splot x**2/2.0+a*y**2/2.0 w l

Then make a contour plot of this function which is equivalent to plotting the phase curves. [explain why.]

Using gnuplot we have the commands:

unset surface
set contour
set size 0.71,1
set cntrparam levels incremental 0.0,10.0,100.0
set view 0,0
unset clabel
unset key
set isosamples 20,20
splot x**2/2.0+a*y**2/2.0 w l

[Note that C=0 at the fixed point (0,0) C>0 for other contours.]

Exercise (Anharmonic oscillator)
(a) Plot the phase portrait for the anharmonic oscillator
d2x/dt2 = -x + a x3.
(b) Plot the stable and unstable manifolds for the linearized system around one of the saddle points.

(a) Set y=dx/dt, thus we have the system of equations
dx/dt=y
dy/dt=-x + a x3

Fixed Points:
(i) x=0, y=0
(ii) x=1/sqrt(a), y=0
(iii) x=-1/sqrt(a), y=0

Phase curves are given by
dy/dx = (-x + a x3)/y ==> y2 + x2/2 - a x4/4 = C

Let us define the function f(x,y)=y2/2 + x2/2 - a x4/4
and plot it. Choose, for example a=0.01, so that fixed points are (0,0), (10,0), (-10,0).
The graph of f(x,y) can be plotted using the gnuplot commands:

a=0.1
set xrange[-15:15]
set yrange[-10:10]
splot y**2/2.+x**2/2.-a*x**4/4. w l

Then, make a contour plot of f(x,y) which is equivalent to plotting the phase curves. [explain why.]
Using gnuplot we have the commands:

unset surface
set contour
set cntrparam levels incremental 0,5,70
set view 0,0
unset clabel
unset key
splot y**2/2.+x**2/2.-a*x**4/4. w l

[Note that C=0 at the fixed point (0,0) while the curve with C=25 contains the points (+-10,0).]

A smoother graph is obtained if we use:

set isosamples 50
splot y**2/2.+x**2/2.-a*x**4/4. w l

Finally, to print your graphs you may use the following commands:

set term post
set output 'phaseportrait.ps'
replot

The graph would look better if we use:

set grid
set xtics -15,5,15
set ytics -10,5,10
set xtics font "Helvetica, 14"
set ytics font "Helvetica, 14"
set title "Simple pendulum - phase portrait" offset 0,-4 font "Helvetica,20"
splot y**2/2.+x**2/2.-a*x**4/4. w l

Here is the result:
Anharmonic oscillator - Phase portrait

(b)

Let us choose a=1. The linearized equations around the fixed point (1,0) are dξ/dt = y, dy/dt = 2ξ, where ξ=x-1. Example for code:

set xrange[0.5:1.5]; set yrange[-0.5:0.5]
set cntrparam levels incremental 0,0.125,0.5
set isosamples 100,100
set title 'Anharmonic oscillator: phase curves (solid), linearized eqns (dashed)' font "Helvetica,16"
splot x**2/2+y**2/2-x**4/4 w l, (x-1)**2/2-y**2/4 w l

Here is the result:
Anharmonic oscillator - Stable and unstable manifolds

Exercise (simple pendulum)
Plot the phase portrait for the simple pendulum
d2x/dt2 + a sin(x) = 0.

Set y=dx/dt, thus we have the system of equations
dx/dt=y
dy/dt=-a sin(x)

Fixed Points:
x=kπ, y=0

Phase curves are given by
(dx/dt)2 - a cos(x) = E, (E:const.)
Using commands similar to Exercise "Anharmonic oscillator" we obtain the result:
Lotka-Voltera solution

Mechanics

Motion in a potential (Lennard-Jones)
We consider the Lennard-Jones potential for chemical bonds
V(r) = D(R/r)12-2(R/r)6.
(a) Plot the potential.
(b) Find the period of oscillations for small deviations from the minimum.
(c) Find the period of oscillations as a function of energy.

(a) The potential has a minimum at r0=R D1/6, where V(r0)=-1/D.
(see lecture notes for further details).

We may use matlab to plot the graph of the potential:

D= 1;
R= 5;
r=0:0.2:10; --> generate a series of points in the r-axis
V=D*(R./r).^12-2*(R./r).^6; --> calculate the corresponding V(r)
plot(r,V) --> plot the graph
axis([0 10 -1 1]) --> define the range of the axes
xlabel('distance') --> label of horizontal axis
ylabel('potential') --> label of vertical axis

(b) We find the second derivative of V at the minimum: k:= d2V(r0)/dr2 = 72/(R2 D4/3). Therefore we may approximate the potential near its minimum by
V(r) = -1/D + (1/2) k (r-r0)2,
The angular frequency of oscillations is ω0=sqrt(k) and the period of oscillations is T0=2π/sqrt(k)=....

To visualize the approximation we may use the following

r0=R*D^(1/6)
V0=-1/D
k=72/(R^2*D^(4/3))
V1 = V0 + 0.5*k*(r-r0).^2;
hold on
plot(r,V1,'--r'); --> plots red dashed line

(c) Let the energy have the value E0 (e.g., E0=-0.5, E0=0.0, etc).
The motion will take place at the r where E0-V(r) >= 0.
If we find the points r1, r2 where E0-V(r) = 0, we then infer that the motion takes place in the interval r1 < r < r2. The solutions of the above equation is found by using the commands:

global D R E0;
E0 = -0.5;
r1 = fzero('len_jon',4) --> solve the equation V-E0=0, using initial guess r=4
r2 = fzero('len_jon',6) --> use different initial guess, thus finds the second solution

The above code uses the function:

function lj=len_jon(r)
global D R E0;
V = D*(R/r)^12-2*(R/r)^6;
lj = E0 - V;

The period of oscillation T is given by the integral:
T/2 = int{ dr/[2 Sqrt(E0-V(r))] } (with limits of integration r=r1, r=r2).
This is calculated using the command:

period = 2*quad('period_lj',r1,r2) --> numerical calculation of the integral of the function period_lj in the interval [r1,r2]

Note: you may need to write: r1=r1+0.0001, r2=r2-0.0001 (for the integral to converge numerically).

The above uses, as an integrand, the function:

function lj=period_lj(r)
global D R E0;
V = D*(R./r).^12-2*(R./r).^6;
lj = 1./sqrt(2*(E0-V));

We can repeat the procedure for other values of Ε0. Use the code

i=0;
for E0=-0.98:0.02:-0.3
i=i+1;
energy(i) = E0;
r1 = fzero('len_jon',4);
r2 = fzero('len_jon',6);
r1 = r1+0.0001;
r2 = r2-0.0001;
period(i) = 2*quad('period_lj',r1,r2);
end

This produces T=T(E0). and we plot the result:

plot(energy,period,'o');
axis([-1 0 0 10]);

Here is the graph:
period versus energy for Lennard-Jones potential