Mathematical Modeling: laboratory exercise 1
(a) graphs and (b) solutions to algebraic equations
Gnuplot
Use gnuplot, a command-line driven graphing utility to
plot graphs, contours, and make surface plots using gnuplot.
Online help is available in, e.g.,
Exercise 1
Suppose a mass m attached to a spring with konstant k.
Then we know that we have harmonic motion with frequency ω=square-root(k/m),
that is, it has a period T=2 π/ω.
Let us plot the function x(t).
m = 1.0
k = 1.0
omega = sqrt(k/m)
period = 2*pi/omega
set xrange[0:period]
set yrange[-1:1]
set xlabel 'time'
set ylabel 'x(t)'
plot 0.9*cos(omega*x)
Γιά να αποθηκεύσουμε το γράφημα μας σε φάκελο postscript με όνομα
fig.ps δίνουμε τις εντολές:
set term post
set output 'fig.ps'
replot
Γιά να δούμε πάλι το γράφημα μας στην οθόνη:
set term x11
replot
Exercise 2
Let us plot the Lennard-Jones potential
V(r) = D(R/r)12-2(R/r)6.
Recall that it has a minimum at r0=D1/6R
with V(r0)=-1/D.
D = 1.0
R = 5.0
r0 = D**(1./6.)*R (we expect the function minimum at r0)
set xrange[0:2*r0]
set yrange[-1.2:1]
set xlabel 'atom separation'
set ylabel 'V(r)'
plot D*(R/x)**12 - 2.0*(R/x)**6
We also know that its second derivative is
d2V(r0)/dr2 = 72/(R2 D4/3).
Therefore we may approximate the potential near its minimum by
V0 = -1/D
k = 72.0/(R**2*D**(4./3.))
replot V0+0.5*k*(x-r0)**2
Matlab
Use matlab to plot graphs and find solutions for algebraic equations and
for the initial value problem of differential equations.
Online help is available in, e.g.,
Exercise 1
A mass m=1 is attached to a rod of length l=1
(the acceleration of gravity is assumed g=1).
The potential of the system is V = -cosθ
(a) Plot the potential.
x=-pi:0.02:pi; --> generate a series of points in the x-axis
y=-cos(x); --> calculate the corresponding y(x)
plot(x,y) --> plot the graph
axis([-pi pi -1 1]) --> define the range of the axes
xlabel('angle') --> label of x-axis
ylabel('potential') --> label of y-axis
(b) Let the energy have the value E0
(e.g., E0=-0.9, E0=0.0)
and find the interval where the motion will take place.
The motion will take place at the x where
E0-V(x) >= 0.
We first need to find the points x1, x2 where
we have E0-V(x) = 0.
The we infer that the motion takes place in the interval
x1 < x < x2.
The solutions of the above equation is found by using the commands:
x1 = fzero('-0.9+cos(x)',-0.1) --> solve the equation -0.9+cos(x)=0, using initial guess x=0.1
x2 = fzero('-0.9+cos(x)',0.1) --> use different initial guess, thus finds the second solution
(c) Find the period of oscillation T.
(Compare this to the period of oscillation for harmonic motion
T=2π).
The period of oscillation is given by the integral:
T/2 = int{ dx/[2 Sqrt(E0+cos(x))] } (with limits of integrtion x=x1, x=x2).
This is calculated using the commands:
f= inline('1./sqrt(2*(-0.9+cos(x)))') --> define the integrand
period = 2*quad(f,x1,x2) --> numerical calculation of the integral
We may repeat the procedure for Ε0=0.0, etc.
We may produce T=T(E0).
(d) Solve the differential equation of motion and plot the solution.
We first have to write a file
pendulum.m
which contains the equatios of motion (i.e., the time derivatives of
the variables).
Then we should tell matlab where this file is located,
e.g., if it is in a directory called newpath then we type the command:
path(path,'newpath')
The differential equation is solved by using the commands:
t0=0;
tfinal=2.0*pi;
y0 = [0 1]'; --> initial conditions
[t,y] = ode45('pendulum',[t0 tfinal],y0); --> numerical solution using Runge-Kutta method
We finally plot the results by:
plot(y(:,1),y(:,2),'-',y(:,3),y(:,4),'-'); --> this needs fixing
Help on the commands used in this exercise can be found in the matlab documentation: