Code covered by the BSD License  

Highlights from
5 ways to Solve IVP ODEs in MATLAB

image thumbnail

5 ways to Solve IVP ODEs in MATLAB

by

 

The scripts show 5 ways of solving IVP ODEs in MATLAB environment

IVP_ODEsols_5_ways.m
%% Solving IVP ODEs using Symbolic MATH "dsolve", Laplace Transforms and
% MuPAD for analytic solutions and ODE45, ODE113 and Simulink for numeric
% solutions of the given 2nd order ODE is shown in this script. 
%
% Given a 2nd order ODE: 3*y''+3*y'+9*y=15*t with ICs: y(0)=3, y'(0)=-2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                by Sulaymon L. ESHKABILOV, Ph.D
%                         December, 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
clear
%% 1st Solution method: DSOLVE
Solution_dsolve=dsolve('3*D2y+3*Dy+9*y-15*t', 'y(0)=3', 'Dy(0)+2=0');
display('Symbolic MATH dsolve function output is: ')
disp(Solution_dsolve)
% display plot of the solution of dsolve and compare it with LT
figure;
subplot(211)
ezplot(Solution_dsolve)
title('Solution found with DSOLVE '); grid
%% 2nd Solution Method: Laplace Transforms 
% Define symbolic variables' names
syms t s Y
ODE2nd='3*D(D(y))(t)+3*D(y)(t)+9*y(t)-15*t';

% Laplace Transforms
lt_A=laplace(ODE2nd, t, s);

% Substitute Initial Conditions and initiate an unknown Y
lt_A=subs(lt_A,{'laplace(y(t),t,s)','y(0)','D(y)(0)'},{Y,3,-2});

% Solve for Y unknown
Y=solve(lt_A, Y);

display('Laplace Transforms of the given 2nd Order ODE with ICs')
disp(Y)

Solution_Laplace=ilaplace(Y);
display('Solution found using Laplace Transforms ')
disp(Solution_Laplace)
subplot(212)
ezplot(Solution_Laplace); grid
title('Solution found with Laplace Transforms')
hold off
 
 % Plot two solutions in one plot area
 
 t=0:pi/100:6*pi;
 SDsol=eval(vectorize(Solution_dsolve));
 LTsol=eval(vectorize(Solution_Laplace));
 figure;
 plot(t, SDsol, 'bx-'); grid
 hold on; xlabel('t'); ylabel('solution values')
 plot(t, LTsol, 'ro-');
 
%% 3rd Solution Method: MuPAD
% MuPAD commands are: 
% % % ooO:=ode({3*y''(t)+3*y'(t)+9*y(t)=15*t, y(0)=3, y'(0)=-2}, y(t))
% % % solutions:=ode::solve(ooO)
% % % plot(solutions, #G, #L)

%% 4th Solution Method: ODE45, ODE113 solvers to get numeric solutions
% Help: to get this cell executed the following function file needs to be
% created in separate (e.g. compareODEsols.m):
% % % % % function dydt=compareODEsols(t, y)
% % % % % % This is a 2nd order ODE IVP example
% % % % % %    3*y''+3*y'+9*y=15*t with ICs: y(0)=3, y'(0)=-2;
% % % % % 
% % % % % dydt=[y(2);-(y(2)+3*y(1))+5*t ];
% % % % % return

ICs=[3, -2];
timeSPAN=0:pi/100:6*pi;
options_ODE=odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
[t, y_ODE45]=ode45(@compareODEsols, timeSPAN, ICs, options_ODE);
plot(t, y_ODE45(:,1), 'm+-');
[t, y_ODE113]=ode113(@compareODEsols, timeSPAN, ICs, options_ODE);
plot(t, y_ODE113(:,1), 'k<-');
title('Comparison analysis of ODE sols with ODE45, ODE23, dsolve, Laplace Transforms & Simulink')


%% 5th Solution Method: Simulink 
% NB: SIMULINK models uses ODE solver ODE23s in simulation of its *.mdl
% model. A solver type can be altered if needed for a better/faster solution search.
% This cell first recalls and executes Simulink model ODEsols_COMPARE.mdl
% that is located and it to be put in currect directory of MATLAB from the archived folder.
sim('ODEsols_COMPARE.mdl');
SIM_sols=ODEsols_SIM.signals.values; time=tout;
plot(time, SIM_sols, 'cd-'); 
legend('dsolve','Laplace Transforms','ODE45','ODE113','Simulink', 0)
hold off
shg

Contact us