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
Let \( F_N (x) \) be the finite Fourier sum with N+1 terms: