Mathematical Modeling: laboratory exercise 2

solution to differential equations: initial value problem

Exercise 1

Suppose two vortices of stength Q1 and Q2 interacting with a potential
V(l) = -Q1 Q2 ln(l)
where l is the distance between the vortices. (a) Solve the equations of motion for Q1, Q2 and for initial conditions of your choice. (b) Plot the coordinates of the vortices as a function of time. (c) Plot the trajectories of the vortices.
(a) Solve the equations of motion.
We use matlab to solve this problem. We first need to write a file called, e.g., vortices.m which should contain the four equations of motion, that is the first time derivatives of the four variables, which are the four coordinates of the vortices (see the lecture notes for the equations).
Note: we should tell matlab where the file containing the function (vortices.m) is located. If this is in a directory called newpath then we should use the command: path(path,'newpath').

We may then use the commands:

t0 = 0; (initial time)
t1 = 20; (final time)
y0 = [-1 0 1 0]; (y0 is the array containing the initial vortex positions, (-1,0) and (1,0) in this case)
[t,y] = ode45('vortices',[t0,t1],y0); (integrates using Runge-Kutta 4th-5th order, from time=t0 to t1, and initial condition y0)

The result of the integration is stored in array y at times which are stored in array t.

Note: compare the result with the result from a more accurate numerical method (or even with a less accurate numerical method).
(b) Plot the coordinates with time.
A plot of the first coordinate as a function of time is obtained by the command:
plot(t(:),y(:,1),'-');

We may also make the plot more clear by using various commands, e.g.,:
axis([t0 t1 -1 1]); (define the range of the axes)
xlabel('time'); (label of x-axis)
ylabel('vortex coordinates'); (label of y-axis)
(c) Plot the trajectories.
A plot of the trajectories is obtained by the command
plot(y(:,1),y(:,2),'-',y(:,3),y(:,4),'-');

Note: set the plot window so that the units of length in both axes (x,y) are equal.


Help on the commands used in this exercise can be found in the matlab documentation:

Exercise 2

Suppose the Lennard-Jones potential for chemical bonds:
V(r) = D(R/r)12 - 2 (R/r)6,
whete r is the distance between two atoms and D, R are constants. Let the equation of motion be d2r/dt2=-dV/dr, then (a) plot the potential, (b) find the period of the motion near the equalibrium point (harmonic approximation), (c) find the period of motion as a function of the energy of the system.
(a) Plot of the potential.
We may use, for example, the commands:
d=1.0;
r=5.0;
x=4:0.1:9;
y = d*r^12*x.^(-12) -2*r^6*x.^(-6);
plot(x,y)

Check that the minimum is at r0=D1/6 R and it has the value V(r0) = -1/D.
(b) Period of motion near equilibrium.
The frequency is ω0 = sqrt(72)/(R D2/3) (see lecture notes), thus the period is T=2π/ω.
(c) Period of motion as a function of energy.
Follow the same procedure as in laboratory exercise 1.
We will now need the function

f = inline('1.0./sqrt(2.0.*(e0-d.*(r./x).^12+2.0.*(r./x).^6))')
where. for example,
d=1.0,
r=5.0, and
e0=-0.5 (e0 is the energy and we should have -1/d < e0 < 0)

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