Prof. Vladimir A. Dobrushkin
This tutorial contains many Matlab scripts.
You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to
I. Verification
function solveDifferentialEquation
%This function solves y'+y = x from starttime to endtime with y(0) = intialCondition
%Parameter initializations
starttime = 0;
endtime = 20;
initialCondition = -2;
function yprime = DifferentialEquationFunction(x,y)
yprime = y - x;
end
[xVals,yVals] = ode45(@DifferentialEquationFunction,[starttime,endtime],initialCondition);
plot(xVals,yVals);
end %Function Complete
function solveDifferentialEquation2
% Function to solve the differential equation
% of a particle in motion due to impact of an external sinusoidal
% force ( yprime = F*sin(w*t) - C*y )
% from a starttime to stoptime, where the intial y value is assigned as
% initialCondition
%Parameter initializations
n = 5;
onesv = [1, 1.45, 1.34, 1.81, 3.5]; %acts as a scalar vector that scales the parameters by its %entries.
C = 10.*onesv;
F = 1.*onesv;
w = 1.*onesv; %With scaled frequencies,w, the plots should show a trend of increasing spread with t5he 3rd plot as an exception since 1.34 defies the pattern.
starttime = 0;
endtime = 20;
initialCondition = 0;
counter = 1;
%Loop to solve the ODE and plot five different graphs.
for indeX = 1:n
[t,y] = ode45(@DifferentialEquationFunction2,[starttime,endtime],initialCondition);
counter = counter + 1;
figure(indeX)
plot(t,y);
title('Solution to differential equation')
xlabel('time')
ylabel('y-value')
end
function yprime = DifferentialEquationFunction2(t,y) % Function defining the ODE
yprime = F(counter)*sin(w(counter)*t) - C(counter)*y;
end
end %Function Complete
Example 1.1.1:
Example 1.1.2:
II. Plotting
The best way to plot direction fileds is to use existing m-files, credited to John Polking from Rice University (http://math.rice.edu/~polking/). He is the author of two special MATLAB routines: dfield8 plots direction fields for single, first order ordinary differential equations, and allows the user to plot solution curves; pplane8 plots vector fields for planar autonomous systems. It allows the user to plot solution curves in the phase plane, and it also enables a variety of time plots of the solution. It will also find equilibrium points and plot separatrices. The latest versions of dfield8 and pplane8 m-functions are not compatible with the latest Matlab version (R2015b). You can grab its modification that works from dfield.m
Since Matlab has no friendly subroutine to plot direction fileds for ODEs, we present several codes that allow to plot these fields directly. Many others can found on the Internet.
function dirfield(f,tval,yval)
% dirfield(f, t1:dt:t2, y1:dy:y2) from
% http://terpconnect.umd.edu/~petersd/246/dirfield.m
%
% plot direction field for first order ODE y' = f(t,y)
% using t-values from t1 to t2 with spacing of dt
% using y-values from y1 to t2 with spacing of dy
%
% f is an @ function, or an inline function,
% or the name of an m-file with quotes.
%
% Example: y' = -y^2 + t
% Show direction field for t in [-1,3], y in [-2,2], use
% spacing of .2 for both t and y:
%
% f = @(t,y) -y^2+t
% dirfield(f, -1:.2:3, -2:.2:2)
[tm,ym]=meshgrid(tval,yval);
dt = tval(2) - tval(1);
dy = yval(2) - yval(1);
fv = vectorize(f);
if isa(f,'function_handle')
fv = eval(fv);
end
yp=feval(fv,tm,ym);
s = 1./max(1/dt,abs(yp)./dy)*0.35;
h = ishold;
quiver(tval,yval,s,s.*yp,0,'.r'); hold on;
quiver(tval,yval,-s,-s.*yp,0,'.r');
if h
hold on
else
hold off
end
axis([tval(1)-dt/2,tval(end)+dt/2,yval(1)-dy/2,yval(end)+dy/2])
Enter text here
function project
% initialising constants
g = 32;
l = 41;
k = 10;
s0 = (pi/2.1); % s0 - initial value of theta
sp0 = 0; % sp0 - initial value of theta prime
trange1 = [0:0.01:13.15]; % time range
w0 = [s0, sp0]; % w0 - initial values
[t, w] = ode45(@solver,trange1,w0);
[t2, w2] = ode45(@solver2,trange1,w0);
%assign theta from both functions
s1 = w(:,1);
s2 = w2(:,1);
%calculate theta double prime for both functions
spp1 = ((-(k^2)*(cos(s1)))./((l^2)*((sin(s1)).^3)));
spp2 = (((-(3*g))/(2*l))*(cos(s2)));
% plot both results
figure
plot(s1,spp1)
hold on
plot(s2,spp2)
legend('contact with wall', 'freefall')
xlabel('\theta')
ylabel('\theta"')
title('Transition of a sliding ladder to a falling ladder')
%set min and max bounds
axis([0 1.5 -1.5 1])
function dwdt=solver(t, w)
s=w(1);
sp=w(2);
dwdt = [sp; ((-(k^2)*(cos(s)))/((l^2)*((sin(s))^3)))];
end
function dwdt=solver2(t, w)
s=w(1);
sp=w(2);
dwdt = [sp; ((-(3*g)/(2*l))*(cos(s)))];
end
end