Code covered by the BSD License  

Highlights from
ODE Tutorial Demo

image thumbnail

ODE Tutorial Demo

by

 

08 Nov 2011 (Updated )

Advanced examples of using the MATLAB Symbolic Math Toolbox for Ordinary Differential Equations

ODEtutorial.m
% This demo shows examples of using the MATLAB Symbolic Math Toolbox for 
% Ordinary Differential Equations (ODE). The examples range from relatively 
% basic to rather advanced, and are placed according to the topics, as
% traditionally covered in undergraduate ODE classes. 
%
% The recommended use of this script is to run it line-by-line from the
% MATLAB editor in the debug mode. Or just copy-paste the relevant part of
% the script into the MATLAB command window.

% License: BSD 
% Copyright (c) 2011 A.V. Knyazev, Andrew.Knyazev@ucdenver.edu
% http://math.ucdenver.edu/~aknyazev/
% $Revision: 1.1 $  $Date: 7-Dec-2011
% Tested in MATLAB 7.13 (R2011b) and its Symbolic Math Toolbox 5.7 (R2011b) 
% Does NOT work in Octave 3.4.2 and below because of poor symbolic toolbox

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Derivatives, differential equations, and initial conditions');
disp(' '); clear all; close all; 
disp('Define symbolic variables, using SYM, and function y=C*exp(t^2)'); 
t=sym('t'); C=sym('C'); y = C*exp(t^2); disp(y); 
disp('The default notation for an independent variable is t.');
disp('Calling diff(y) takes the first derivative of y with respect to t,');
yprime = diff(y); disp(yprime); 
disp('so the result is exactly the same as diff(y,t)');
yprime = diff(y,t); disp(yprime); 
disp('Calling diff(y,2) takes the second derivative of y with respect to t,');
y2prime = diff(y,2); disp(y2prime); 
disp('Directly check that y=C*exp(t^2) solves D2y-2*t*Dy-2*y=0');
if y2prime-2*t*yprime-2*y == 0
    disp('Checked!');
else
    disp('Something is wrong!');
end
disp('Take some value for C, e.g., C=2, substitute in y, using SUBS,');
C=2; yC2=subs(y); disp(yC2); 
disp('and plot the resulting function.');
ezplot(yC2,[-1,1]);
disp('Finally, pick up some value for t, e.g., t=1, and evaluate y, using SUBS,');
yC2t1=subs(yC2,'t',1); display(yC2t1); % which is the same as 
t=1; yC2t1=subs(yC2); display(yC2t1); 
disp('Press any key to continue to the next section'); 
pause; clear all; close all;  disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Plotting <a href = "http://en.wikipedia.org/wiki/Slope_field">Slope Field</a>');
disp('Press any key to continue'); 
pause; clear all; close all; 
disp('Plot the slope field and isoclines for f(x,y)=xy');
Ffun = @(X,Y)X*Y;                % in-line function f(x,y)=xy
Ffun=str2func(vectorize(Ffun));  % inserts a '.' before '*' for vectors
                                 % post 2011b MATLAB may not need str2func
[X,Y]=meshgrid(-2:.3:2,-2:.3:2); % choose the plot sizes
DY=Ffun(X,Y); DX=ones(size(DY)); % generate the plot values 
quiver(X,Y,DX,DY);               % plot the direction field
hold on;
contour(X,Y,DY,[-6 -2 -1 0 1 2 6]); %add the isoclines
title('Slope field and isoclines for f(x,y)=xy');
% the code above also works in GNU Octave
disp('Press any key to continue to the next section'); 
pause; clear all; close all; disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
link = 'http://en.wikipedia.org/wiki/Autonomous_differential_equation';
linktitle = 'Autonomous Differential Equations';
fprintf('First-order <a href = "%s">%s</a>\n', link, linktitle) 
disp('The equation Dy = (2-y)y is autonomous, since the independent variable,'); 
disp('let us call it x, does not explicitly appear in the equation.');
disp('To plot the slope field and isocline for this equation, ');
disp('one can use the following code in GNU Octave/MATLAB.');
disp('Press any key to continue'); pause; clear all; close all; 
Ffun = @(X,Y)(2-Y)*Y;            % in-line function f(x,y)=(2-y)y
Ffun=str2func(vectorize(Ffun));  % inserts a '.' before '*' for vectors
[X,Y]=meshgrid(0:.2:6,-1:.2:3);  % choose the plot sizes
DY=Ffun(X,Y); DX=ones(size(DY)); % generate the plot values 
quiver(X,Y,DX,DY);               % plot the direction field
hold on;
contour(X,Y,DY,[0 2]); %add the isoclines
title('Slope field and isoclines for f(x,y)=(2-y)y')
disp('One can observe from the plot that the function (2 - y)y is, '); 
disp('of course, x-invariant, and so it the shape of the solution, i.e.,'); 
disp('y(x) = y(x-x0) for any shift x0.');
disp(' ');
disp('Now we solve the equation symbolically in MATLAB, by DSOLVE');
y=dsolve('Dy=(2-y)*y','x'); % solve the equation symbolically
disp(y); %shows the solution y 
disp('MATLAB finds two equilibrium solutions, y = 0 and y = 2,'); 
disp('and a third solution involving an unknown constant C3,');
disp('y(3)=-2/(exp(C3 - 2*x) - 1)');
disp(' ');
disp('Piking up some specific values for the initial condition,'); 
disp('we can add the plot of several solutions.');
y1=dsolve('Dy=(2-y)*y','y(1)=1','x'); % solve the initial value problem symbolically 
y2=dsolve('Dy=(2-y)*y','y(2)=1','x'); % for different initial conditions 
y3=dsolve('Dy=(2-y)*y','y(3)=1','x'); y4=dsolve('Dy=(2-y)*y','y(1)=3','x'); 
y5=dsolve('Dy=(2-y)*y','y(2)=3','x'); y6=dsolve('Dy=(2-y)*y','y(3)=3','x');
ezplot(y1, [0 6]); ezplot(y2, [0 6]); % plot the solutions
ezplot(y3, [0 6]); ezplot(y4, [0 6]); ezplot(y5, [0 6]); ezplot(y6, [0 6]); 
title('Slope field, isoclines and solutions for f(x,y)=(2-y)y')
disp('Press any key to continue to the next section'); 
pause; clear all; close all; disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
link = 'http://en.wikipedia.org/wiki/Ordinary_differential_equation#Linear_ordinary_differential_equations';
linktitle = 'First-order linear homogeneous ordinary differential equations';
fprintf('First-order <a href = "%s">%s</a>\n', link, linktitle) 
disp('Press any key to continue'); pause; clear all; close all; 
disp('Solve equation Dy+f(t)y=0 symbolically in MATLAB, by DSOLVE');
y = dsolve('Dy + f(t)*y=0'); % solve the equation symbolically
disp(y); %shows the solution y = C1/exp(int(f(t), t))
disp('Replace function f(y) with a particular choice f(t)=sin(t) in the solution');
disp('can be done by the SUBSEX function, but not by SUBS');
ysin = feval(symengine, 'subsex', y, ['f(t)=sin(t)']); 
ysin=simplify(ysin); disp(ysin);
disp(' '); disp('Press any key to continue'); pause; clear all; close all;
disp('Solve equation Dy+ay=0, and find constant a, so y(1)=2 and y(2)=1');
disp('By default, MATLAB symbolic toolbox calls the independent variable "t".');
disp('There is no even need to declare t as symbolic using SYM. We get'); 
% t=sym('t');  % declare a symolic variable t - not necessary
y = dsolve('Dy + a*y=0', 'y(1)=2', 'y(2)=1'); % solve the equation symbolically
disp(y); %shows the solution y = (2*exp(a))/exp(a*t); 
disp('Now, solving symbolically the nonlinear equation for "a" by SOLVE:');
a=solve('2*exp(a) = exp(2*a)'); disp(a);
disp('we symbolically plug the just determined value for "a" into "y" by SUBS:')
z=subs(y); disp(z); % we obtain the solution 4/exp(t*log(2))
disp('But this explicit solution can still be simplified by SIMPLIFY')
simplify(z)
disp(' '); disp('Press any key to continue'); pause; clear all; close all;
disp('Solve the initial value problem t*Dy+e*y=0, y(-e)=e, where e=exp(1)'); 
disp('Note that MATLAB uses the LOG function for Natural logarithm.'); 
disp('We solve by DSOLVE and simplified by SIMPLIFY:')
y=dsolve('t*Dy+e*y=0','y(-e)=e'); y=simplify(y); disp(y); % -(-e)^(e + 1)/t^e
disp('In order to check the initial condition and plot the solution,');
disp('the constant "e" needs to be defined and substituted.')
e=sym('exp(1)'); y=subs(y);
disp('Let us check the initial condition:'); 
t=-e; y0=subs(y); checkBC = simplify(y0)-e; 
if checkBC == 0
    disp('The initial condition is satisfied');
else
    disp('The initial condition is NOT satisfied. Something is wrong!');
end%if
disp('Plotting the solution by EZPLOT'); ezplot(y);
disp('Press any key to continue to the next section'); 
pause; clear all; close all; disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
link = 'http://en.wikipedia.org/wiki/Ordinary_differential_equation#Linear_ordinary_differential_equations';
linktitle = 'First-order linear non-homogeneous ordinary differential equations';
fprintf('First-order <a href = "%s">%s</a>\n', link, linktitle)
disp('Press any key to continue'); pause; clear all; close all; 
disp('Solve equation Dy+p(t)y=g(t) symbolically in MATLAB, by DSOLVE');
y = dsolve('Dy + p(t)*y=g(t)'); disp(y); % results in
% y = C1/exp(int(p(t), t)) + int(exp(int(p(t), t))*g(t), t)/exp(int(p(t), t))
disp('Press any key to continue to the next section'); 
pause; clear all; close all; disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
link = 'http://en.wikipedia.org/wiki/Examples_of_differential_equations#Separable_first-order_ordinary_differential_equations';
linktitle = 'First-order nonlinear separable ordinary differential equations';
fprintf('First-order <a href = "%s">%s</a>\n', link, linktitle)
disp('Press any key to continue'); pause; clear all; close all; 
disp('Solve equation (2*y-sin(y))*Dy+t=sin(t), y(0)=0 symbolically in MATLAB, by DSOLVE');
y=dsolve('(2*y-sin(y))*Dy+t=sin(t)', 'y(0)=0'); disp(y); % results in
% solve(- cos(y) - y^2 = cos(t) + t^2/2 - 2, y) % implicit solution
link = 'http://www.wolframalpha.com/input/?i=%282y-sin%28y%29%29y%27%2Bt%3Dsin%28t%29%2C+y%280%29%3D0';
linktitle = 'solution by WolframAlpha';
fprintf('Compare to <a href = "%s">%s</a>\n', link, linktitle);
disp('Press any key to continue to the next section'); 
pause; clear all; close all; disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
linktitle = 'Bernoulli differential equation'; 
link = 'http://en.wikipedia.org/wiki/Bernoulli_differential_equation'; 
fprintf('First-order <a href = "%s">%s</a>\n', link, linktitle)
disp('Press any key to continue'); pause; clear all; close all; 
disp('Solve equation Dy-2*y/x=-x^2*y^2 symbolically in MATLAB, by DSOLVE');
x = dsolve('Dy-2*y/x=-x^2*y^2','x') ; disp(x); 
% gives both solutions: 0 and x^2/(x^5/5 + C1)
linktitle = 'solution by WolframAlpha';
link = 'http://www.wolframalpha.com/input/?i=y%27-2*y%2Fx%3D-x%5E2*y%5E2';
fprintf('Compare to <a href = "%s">%s</a>\n', link, linktitle)
disp('Solve equation Dy=2*t*y*(1-y) symbolically in MATLAB, by DSOLVE');
y=simplify(dsolve('Dy=2*t*y*(1-y)', 'y(0)=-1')); disp(y);
%gives -1/(2/exp(t^2) - 1)
disp('and determine the points where the denominator is zero, by SOLVE');
singularities=solve('2/exp(t^2)=1'); disp(singularities); 
% log(2)^(1/2) and - log(2)^(1/2)
linktitle = 'solution by WolframAlpha';
link = 'http://www.wolframalpha.com/input/?i=y%27%3D2*t*y*%281-y%29%2C+y%280%29%3D-1';
fprintf('Compare to <a href = "%s">%s</a>\n', link, linktitle);
disp('Press any key to continue to the next section'); 
pause; clear all; close all; disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Second-order linear ordinary differential equations');
disp('Press any key to continue');
pause; clear all; close all; 
linktitle = 'solution by WolframAlpha';
link = 'http://www.wolframalpha.com/input/?i=x%27%27%3D-x';
fprintf('Solve D2x=-x and compare to <a href = "%s">%s</a>\n', link, linktitle);
x = dsolve('D2x=-x','t'); disp(x); % gives x = C1*cos(t) + C2*sin(t)
disp('Press any key to continue'); pause; clear all; close all; 
linktitle = 'solution by WolframAlpha';
link = 'http://www.wolframalpha.com/input/?i=x%27%27%3D-x%2Cx%280%29%3D1%2Cx%27%280%29%3D0';
fprintf('SOLVE IVP D2x=-x, x(0)=1, Dx(0)=0 and compare to <a href = "%s">%s</a>\n', link, linktitle);
x = dsolve('D2x=-x','x(0)=1','Dx(0)=0'); disp(x); %gives x = cos(t)
disp('Press any key to continue'); pause; clear all; close all; 
disp('Calculate the Wronskian for y1=1 and y2=exp(t) at the point t=0.');
t=sym('t'); y1=sym('y1'); y1=1; y2=exp(t);
Wt=det([y1 y2; diff(y1,t) diff(y2,t)]); disp(Wt); % gives exp(t)
Wt0=subs(Wt,t,0); disp(Wt0); % gives 1
disp('Press any key to continue'); pause; clear all; close all; 
disp('Suppose that the functions y1=1 and y2=exp(t) form the fundamental set.'); 
disp('Find the solution satisfying the initial conditions y(0)=2 and Dy(0)=1.'); 
t=sym('t'); y1=sym('y1'); y1=1; y2=exp(t);
Wt=[y1 y2; diff(y1,t) diff(y2,t)]; disp(Wt);   % gives the Wronskian 
Wt0=subs(Wt,t,0); disp(Wt0); % gives the Wronskian at t=0
C = Wt0\[2; 1];   %gives [1 ;1] - the coefficients of the linear combination
Answer = [y1 y2]*C; % gives 1+exp(t)
disp('Press any key to continue'); pause; clear all; close all; 
linktitle = 'Free spring-mass vibrations with dampling';
link = 'http://en.wikipedia.org/w/index.php?title=Examples_of_differential_equations#A_more_complicated_model';
fprintf('<a href = "%s">%s</a>\n', link, linktitle);
disp('Solve IVP D2x+c*Dx+k*x=0, x(0)=1, Dx(0)=0');
x = dsolve('D2x+c*Dx+k*x=0','x(0)=1','Dx(0)=0'); disp(x);
a=sym('a'); b=sym('b'); cs=-2*a; ks=a.^2+b.^2;
a=-0.25; b=1;  c=subs(cs); k=subs(ks); x1=subs(x); 
a=-0.5;  b=1;  c=subs(cs); k=subs(ks); x2=subs(x); 
a=-0.5;  b=.5; c=subs(cs); k=subs(ks); x3=subs(x); 
h3 = ezplot(x3, [0 10]); set(h3, 'Color', 'r'); hold on;
h2 = ezplot(x2, [0 10]); set(h2, 'Color', 'g'); grid on;
h1 = ezplot(x1, [0 10]); set(h1, 'Color', 'b'); title('damped oscillation');
legend([h1 h2 h3],'a=-0.25; b=1;','a=-0.5;  b=1;','a=-0.5;  b=.5;')
disp('Press any key to continue'); pause; clear all; close all; 
disp('Solve the ODE L(y)=D2y-2*t*Dy-2*y=0 by using DSOLVE:');
ysol=simplify(dsolve('D2y-2*t*Dy-2*y=0')); disp(ysol); 
disp('Notice that y=exp(t^2) is a solution in the fundamantal set, since L(y1) is');
t=sym('t');  y1=exp(t^2); rhs=diff(y1,2)-2*t*diff(y1)-2*y1; disp(rhs);
disp('The other solution in the fundamental set is');
disp('function y2=exp(t^2)*erf(t). Let us plug it into L(y2)');
y2=exp(t^2)*erf(t); rhs=diff(y2,2)-2*t*diff(y2)-2*y2; disp(rhs);
disp('The latter does not look like zero, but use SIMPLIFY to check that it is');
rhs=simplify(rhs); disp(rhs);
disp('Satisfying the initial conditions y(1)=1 and Dy(1)=3');
ysolIVP=dsolve('D2y-2*t*Dy-2*y=0','y(1)=1','Dy(1)=3'); disp(ysolIVP);
disp('and plot the resulting function.');
ezplot(ysolIVP,[-1,1]);
%title(['$' latex(ysolIVP) '$'],'interpreter','latex') %title in LATEX
disp('Press any key to continue'); pause; close all; clear all;
disp('Finally, do you like movies? Let us make a movie!');
disp('Satisfy the initial conditions y(1)=1 and Dy(1)=A');
disp('where A varies from 2 to 3, and plot the changing solutions');
ysolIVP=dsolve('D2y-2*t*Dy-2*y=0','y(1)=1','Dy(1)=A');
hf=figure; axis tight; set(gca,'NextPlot','replacechildren');
for j=1:101
    A=2+(j-1)*0.01; yA=subs(ysolIVP);
    ezplot(yA,[-1,1]); 
    title('Solutions of D2y-2*t*Dy-2*y=0, y(1)=1, Dy(1)=A'); drawnow; 
    M(j) = getframe(gcf); %captures entire interior of the current figure 
end
disp('Now you have a movie in MATLAB format that can be played in MATLAB');
disp('Press any key to play the MATLAB movie'); pause; 
hf = figure; axis off; movie(hf,M); 
disp('But wait, you want to send the movie to a friend, who has no MATLAB?');
disp('Let us use  MOVIE2AVI function for the conversion');
movie2avi(M,'MyCoolMovieODE.avi');  
disp('With a file browser in your current directory'); 
disp(pwd) % shows the location of the current directory
disp('locate the file called MyCoolMovieODE.avi'); 
disp('It will play in most video players.');
disp(' ');
disp('It is also possible to produce the AVI-format movie directly.');
disp('Press any key and be patient until your see "DONE"'); 
pause; close all; clear all;
ysolIVP=dsolve('D2y-2*t*Dy-2*y=0','y(1)=1','Dy(1)=A');
aviobj=avifile('MyCoolMovieODEDirect.avi','fps',5,'compression','None'); 
hf = figure('visible','off'); %turns visibility of figure off
for j=1:101
    A=2+(j-1)*0.01; yA=subs(ysolIVP); ezplot(yA,[-1,1]); 
    title('Solutions of D2y-2*t*Dy-2*y=0, y(1)=1, Dy(1)=A');
    aviobj=addframe(aviobj,hf); %adds frames to the AVI file
end
aviobj=close(aviobj); %closes the AVI file
close(hf); %closes the handle to invisible figure
disp('DONE! With a file browser in your current directory'); 
disp(pwd) % shows the location of the current directory
disp('locate the file called MyCoolMovieODEDirect.avi'); 
disp('Press any key to continue to the next section'); 
pause; clear all; close all; disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Matrix calculus');
disp('Press any key to continue');
pause; clear all; close all; 
t=sym('t'); A=[1 2*t; exp(t) sin(t)]; B=[1; exp(t)]; diff(A*B)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Linear systems of ODEs');
disp('Press any key to continue');
pause; clear all; close all; 
linktitle = 'solution by WolframAlpha';
link = 'http://www.wolframalpha.com/input/?i=y%27+%3D+{{1%2C2}%2C{2%2C-2}}.y%2B+{t%2C+sin%28t%29}';
fprintf('Solve Dy1= y1+2*y2+t; Dy2=2*y1-2*y2+sin(t) and compare to <a href = "%s">%s</a>\n', link, linktitle);
y = dsolve('Dy1= y1+2*y2+t; Dy2=2*y1-2*y2+sin(t)'); disp(y); 
disp(y.y1); disp(y.y2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Solving a system of differential equations in matrix form'); 
disp('Y′ = AY + B, where A, B, and Y represent matrices');
disp('Press any key to continue');
pause; clear all; close all;
syms t x y;
A = [1 2; -1 1];
B = [1; t];
Y = [x; y];
sys = A*Y + B; disp(sys); 
%sys = x + 2*y + 1 ; t - x + y
disp('The dsolve function does not accept matrices.'); 
disp('We extract the components of the matrix and convert them to strings');
eq1 = char(sys(1)); disp(eq1);
eq2 = char(sys(2)); disp(eq2);
disp('Use the strcat function to concatenate the left and right sides of the equations'); 
disp('and call the DSOLVE function to solve the system');
[x, y] = dsolve(strcat('Dx = ',eq1), strcat('Dy = ',eq2)); 
disp(x); disp(y);

Contact us