Code covered by the BSD License

# ODE Tutorial Demo

### Andrew Knyazev (view profile)

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
%
% 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.

% 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(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
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(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
linktitle = 'First-order linear homogeneous ordinary differential equations';
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(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
linktitle = 'First-order linear non-homogeneous ordinary differential equations';
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(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
linktitle = 'First-order nonlinear separable ordinary differential equations';
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
disp('Press any key to continue to the next section');
pause; clear all; close all; disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)
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)
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;
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;
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';
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');
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;
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);```