MATLAB lessons N.2 - Enrico Bertolazzi

:open:

function Y = ode1(odefun,tspan,y0,varargin)
%ODE1  Solve differential equations with a non-adaptive method of order 1.
%   Y = ODE1(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates 
%   the system of differential equations y' = f(t,y) by stepping from T0 to 
%   T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
%   The vector Y0 is the initial conditions at T0. Each row in the solution 
%   array Y corresponds to a time specified in TSPAN.
%
%   Y = ODE1(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters 
%   P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...). 
%
%   This is a non-adaptive solver. The step sequence is determined by TSPAN.
%   The solver implements the forward Euler method of order 1.   
%
%   Example 
%         tspan = 0:0.1:20;
%         y = ode1(@vdp1,tspan,[2 0]);  
%         plot(tspan,y(:,1));
%     solves the system y' = vdp1(t,y) with a constant step size of 0.1, 
%     and plots the first component of the solution.   
%

if ~isnumeric(tspan)
  error('TSPAN should be a vector of integration steps.');
end

if ~isnumeric(y0)
  error('Y0 should be a vector of initial conditions.');
end

h = diff(tspan);
if any(sign(h(1))*h <= 0)
  error('Entries of TSPAN are not in order.') 
end  

try
  f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
  msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
  error(msg);  
end  

y0 = y0(:);   % Make a column vector.
if ~isequal(size(y0),size(f0))
  error('Inconsistent sizes of Y0 and f(t0,y0).');
end  

neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);

Y(:,1) = y0;
for i = 1:N-1 
  Y(:,i+1) = Y(:,i) + h(i)*feval(odefun,tspan(i),Y(:,i),varargin{:});
end
Y = Y.';

:open:

%
% Example of ODE solver of MATLAB
%

% set some option for the ODE numerical solver
options = odeset('RelTol',1e-7,'AbsTol',[1e-7 1e-7 1e-7]);

% se initial condition for the problem
Y0 = [ -7.5 -3.6 30 ] ;

% Set a vector specifying the interval of integration, [t0,tf].
% The solver imposes the initial conditions at TSPAN(1),
% and integrates from TSPAN(1) to TSPAN(end). 
% To obtain solutions at specific times (all increasing or all decreasing), 
% use TSPAN = [t0,t1,...,tf].
TSPAN = [ 0:0.01:10 ] ;

% load the "pointer" of the function to FUN
FUN = @lorenz ;

% Runge-Kutta embedded of order 2 and 3 by Bogaci and Shampine
[T1,Y1] = ode23 ( FUN, TSPAN, Y0, options ) ;

% Runge-Kutta embedded of order 4 and 5 by Dormand and Prince (first choice scheme)
[T2,Y2] = ode45 ( FUN, TSPAN, Y0, options ) ;

% Adams-bashforth-Moulton PECE
[T3,Y3] = ode113 ( FUN, TSPAN, Y0, options ) ;

% Numerical Differentiation Formula, similar to Backward Differetiation Formula
% by Gear. Useful for stiff problem.
[T4,Y4] = ode15s ( FUN, TSPAN, Y0, options ) ;

% Rosenbrock of order 2. For stiff problem.
[T5,Y5] = ode23s ( FUN, TSPAN, Y0, options ) ;

% Trapeziodal rule for moderately stiff problem
[T6,Y6] = ode23t ( FUN, TSPAN, Y0, options ) ;

% TR-BDF2 an implit Runge-Kutta where the first stage is the 
% trapezoidal rule and teh second stage is a BDF formula of order 2
[T7,Y7] = ode23tb ( FUN, TSPAN, Y0, options ) ;

%
% plot the functions approximated by various solver
%
plot( T1, Y1(:,1), 'r-',  ...
      T2, Y2(:,1), 'g--', ...
      T3, Y3(:,1), 'b:',  ...
      T4, Y4(:,1), 'c-.', ...
      T5, Y5(:,1), 'm--',  ...
      T6, Y6(:,1), 'y:',  ...
      T7, Y7(:,1), 'k-.' )

:open:

%
% Check Van Der Pol equation
% this is a stiff test
%

% set some option for the ODE numerical solver
options = odeset('RelTol',1e-4) ; %,'AbsTol',[1e-7 1e-7 1e-7]);

% se initial condition for the problem
Y0 = [ 2 0 ] ;

% Set a vector specifying the interval of integration, [t0,tf].
% The solver imposes the initial conditions at TSPAN(1),
% and integrates from TSPAN(1) to TSPAN(end). 
% To obtain solutions at specific times (all increasing or all decreasing), 
% use TSPAN = [t0,t1,...,tf].
TSPAN = [ 0:1:3000 ] ;

% load the "pointer" of the function to FUN
FUN = @vanderpol ;

% This is a stiff problem, the first three solver do not work! 
%[T1,Y1] = ode23   ( FUN, TSPAN, Y0 ) ;
%[T2,Y2] = ode45   ( FUN, TSPAN, Y0 ) ;
%[T3,Y3] = ode113  ( FUN, TSPAN, Y0 ) ;

% those solvers can solve a stiff problem
[T4,Y4] = ode15s  ( FUN, TSPAN, Y0 ) ;
[T5,Y5] = ode23s  ( FUN, TSPAN, Y0 ) ;
[T6,Y6] = ode23t  ( FUN, TSPAN, Y0 ) ;
[T7,Y7] = ode23tb ( FUN, TSPAN, Y0 ) ;

%
% plot the functions approximated by various solver
%
% first component
plot( T4, Y4(:,1), 'c-.', ...
      T5, Y5(:,1), 'm--',  ...
      T6, Y6(:,1), 'y:',  ...
      T7, Y7(:,1), 'k-.' )
      
% create a new window to plot
figure()

% second component
plot( T4, Y4(:,2), 'c-.', ...
      T5, Y5(:,2), 'm--',  ...
      T6, Y6(:,2), 'y:',  ...
      T7, Y7(:,2), 'k-.' )

:open:

function Y = ode1(odefun,tspan,y0,varargin)
%ODE1  Solve differential equations with a non-adaptive method of order 1.
%   Y = ODE1(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates 
%   the system of differential equations y' = f(t,y) by stepping from T0 to 
%   T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
%   The vector Y0 is the initial conditions at T0. Each row in the solution 
%   array Y corresponds to a time specified in TSPAN.
%
%   Y = ODE1(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters 
%   P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...). 
%
%   This is a non-adaptive solver. The step sequence is determined by TSPAN.
%   The solver implements the forward Euler method of order 1.   
%
%   Example 
%         tspan = 0:0.1:20;
%         y = ode1(@vdp1,tspan,[2 0]);  
%         plot(tspan,y(:,1));
%     solves the system y' = vdp1(t,y) with a constant step size of 0.1, 
%     and plots the first component of the solution.   
%

if ~isnumeric(tspan)
  error('TSPAN should be a vector of integration steps.');
end

if ~isnumeric(y0)
  error('Y0 should be a vector of initial conditions.');
end

h = diff(tspan);
if any(sign(h(1))*h <= 0)
  error('Entries of TSPAN are not in order.') 
end  

try
  f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
  msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
  error(msg);  
end  

y0 = y0(:);   % Make a column vector.
if ~isequal(size(y0),size(f0))
  error('Inconsistent sizes of Y0 and f(t0,y0).');
end  

neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);

Y(:,1) = y0;
for i = 1:N-1 
  Y(:,i+1) = Y(:,i) + h(i)*feval(odefun,tspan(i),Y(:,i),varargin{:});
end
Y = Y.';

:open:

function Y = ode2(odefun,tspan,y0,varargin)
%ODE2  Solve differential equations with a non-adaptive method of order 2.
%   Y = ODE2(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates 
%   the system of differential equations y' = f(t,y) by stepping from T0 to 
%   T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
%   The vector Y0 is the initial conditions at T0. Each row in the solution 
%   array Y corresponds to a time specified in TSPAN.
%
%   Y = ODE2(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters 
%   P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...). 
%
%   This is a non-adaptive solver. The step sequence is determined by TSPAN
%   but the derivative function ODEFUN is evaluated multiple times per step.
%   The solver implements the improved Euler (Heun's) method of order 2.   
%
%   Example 
%         tspan = 0:0.1:20;
%         y = ode2(@vdp1,tspan,[2 0]);  
%         plot(tspan,y(:,1));
%     solves the system y' = vdp1(t,y) with a constant step size of 0.1, 
%     and plots the first component of the solution.   
%

if ~isnumeric(tspan)
  error('TSPAN should be a vector of integration steps.');
end

if ~isnumeric(y0)
  error('Y0 should be a vector of initial conditions.');
end

h = diff(tspan);
if any(sign(h(1))*h <= 0)
  error('Entries of TSPAN are not in order.') 
end  

try
  f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
  msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
  error(msg);  
end  

y0 = y0(:);   % Make a column vector.
if ~isequal(size(y0),size(f0))
  error('Inconsistent sizes of Y0 and f(t0,y0).');
end  

neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);
F = zeros(neq,2);

Y(:,1) = y0;
for i = 2:N
  ti = tspan(i-1);
  hi = h(i-1);
  yi = Y(:,i-1);
  F(:,1) = feval(odefun,ti,yi,varargin{:});
  F(:,2) = feval(odefun,ti+hi,yi+hi*F(:,1),varargin{:});
  Y(:,i) = yi + (hi/2)*(F(:,1) + F(:,2));
end
Y = Y.';

:open:

function Y = ode3(odefun,tspan,y0,varargin)
%ODE3  Solve differential equations with a non-adaptive method of order 3.
%   Y = ODE3(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates 
%   the system of differential equations y' = f(t,y) by stepping from T0 to 
%   T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
%   The vector Y0 is the initial conditions at T0. Each row in the solution 
%   array Y corresponds to a time specified in TSPAN.
%
%   Y = ODE3(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters 
%   P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...). 
%
%   This is a non-adaptive solver. The step sequence is determined by TSPAN
%   but the derivative function ODEFUN is evaluated multiple times per step.
%   The solver implements the Bogacki-Shampine Runge-Kutta method of order 3.
%
%   Example 
%         tspan = 0:0.1:20;
%         y = ode3(@vdp1,tspan,[2 0]);  
%         plot(tspan,y(:,1));
%     solves the system y' = vdp1(t,y) with a constant step size of 0.1, 
%     and plots the first component of the solution.   
%

if ~isnumeric(tspan)
  error('TSPAN should be a vector of integration steps.');
end

if ~isnumeric(y0)
  error('Y0 should be a vector of initial conditions.');
end

h = diff(tspan);
if any(sign(h(1))*h <= 0)
  error('Entries of TSPAN are not in order.') 
end  

try
  f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
  msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
  error(msg);  
end  

y0 = y0(:);   % Make a column vector.
if ~isequal(size(y0),size(f0))
  error('Inconsistent sizes of Y0 and f(t0,y0).');
end  

neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);
F = zeros(neq,3);

Y(:,1) = y0;
for i = 2:N
  ti = tspan(i-1);
  hi = h(i-1);
  yi = Y(:,i-1);
  F(:,1) = feval(odefun,ti,yi,varargin{:});
  F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});
  F(:,3) = feval(odefun,ti+0.75*hi,yi+0.75*hi*F(:,2),varargin{:});  
  Y(:,i) = yi + (hi/9)*(2*F(:,1) + 3*F(:,2) + 4*F(:,3));
end
Y = Y.';

:open:

function Y = ode4(odefun,tspan,y0,varargin)
%ODE4  Solve differential equations with a non-adaptive method of order 4.
%   Y = ODE4(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates 
%   the system of differential equations y' = f(t,y) by stepping from T0 to 
%   T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
%   The vector Y0 is the initial conditions at T0. Each row in the solution 
%   array Y corresponds to a time specified in TSPAN.
%
%   Y = ODE4(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters 
%   P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...). 
%
%   This is a non-adaptive solver. The step sequence is determined by TSPAN
%   but the derivative function ODEFUN is evaluated multiple times per step.
%   The solver implements the classical Runge-Kutta method of order 4.   
%
%   Example 
%         tspan = 0:0.1:20;
%         y = ode4(@vdp1,tspan,[2 0]);  
%         plot(tspan,y(:,1));
%     solves the system y' = vdp1(t,y) with a constant step size of 0.1, 
%     and plots the first component of the solution.   
%

if ~isnumeric(tspan)
  error('TSPAN should be a vector of integration steps.');
end

if ~isnumeric(y0)
  error('Y0 should be a vector of initial conditions.');
end

h = diff(tspan);
if any(sign(h(1))*h <= 0)
  error('Entries of TSPAN are not in order.') 
end  

try
  f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
  msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
  error(msg);  
end  

y0 = y0(:);   % Make a column vector.
if ~isequal(size(y0),size(f0))
  error('Inconsistent sizes of Y0 and f(t0,y0).');
end  

neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);
F = zeros(neq,4);

Y(:,1) = y0;
for i = 2:N
  ti = tspan(i-1);
  hi = h(i-1);
  yi = Y(:,i-1);
  F(:,1) = feval(odefun,ti,yi,varargin{:});
  F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});
  F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:});  
  F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});
  Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));
end
Y = Y.';

:open:

function Y = ode5(odefun,tspan,y0,varargin)
%ODE5  Solve differential equations with a non-adaptive method of order 5.
%   Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates 
%   the system of differential equations y' = f(t,y) by stepping from T0 to 
%   T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
%   The vector Y0 is the initial conditions at T0. Each row in the solution 
%   array Y corresponds to a time specified in TSPAN.
%
%   Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters 
%   P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...). 
%
%   This is a non-adaptive solver. The step sequence is determined by TSPAN
%   but the derivative function ODEFUN is evaluated multiple times per step.
%   The solver implements the Dormand-Prince method of order 5 in a general 
%   framework of explicit Runge-Kutta methods.
%
%   Example 
%         tspan = 0:0.1:20;
%         y = ode5(@vdp1,tspan,[2 0]);  
%         plot(tspan,y(:,1));
%     solves the system y' = vdp1(t,y) with a constant step size of 0.1, 
%     and plots the first component of the solution.   

if ~isnumeric(tspan)
  error('TSPAN should be a vector of integration steps.');
end

if ~isnumeric(y0)
  error('Y0 should be a vector of initial conditions.');
end

h = diff(tspan);
if any(sign(h(1))*h <= 0)
  error('Entries of TSPAN are not in order.') 
end  

try
  f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
  msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
  error(msg);  
end  

y0 = y0(:);   % Make a column vector.
if ~isequal(size(y0),size(f0))
  error('Inconsistent sizes of Y0 and f(t0,y0).');
end  

neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);

% Method coefficients -- Butcher's tableau
%  
%   C | A
%   --+---
%     | B

C = [1/5; 3/10; 4/5; 8/9; 1];

A = [ 1/5,          0,           0,            0,         0
      3/40,         9/40,        0,            0,         0
      44/45        -56/15,       32/9,         0,         0
      19372/6561,  -25360/2187,  64448/6561,  -212/729,   0
      9017/3168,   -355/33,      46732/5247,   49/176,   -5103/18656];

B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];

% More convenient storage
A = A.'; 
B = B(:);      

nstages = length(B);
F = zeros(neq,nstages);

Y(:,1) = y0;
for i = 2:N
  ti = tspan(i-1);
  hi = h(i-1);
  yi = Y(:,i-1);  
  
  % General explicit Runge-Kutta framework
  F(:,1) = feval(odefun,ti,yi,varargin{:});  
  for stage = 2:nstages
    tstage = ti + C(stage-1)*hi;
    ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));
    F(:,stage) = feval(odefun,tstage,ystage,varargin{:});
  end  
  Y(:,i) = yi + F*(hi*B);
  
end
Y = Y.';

:open:

function ode_test
%ODE_TEST  Run non-adaptive ODE solvers of different orders.
%   ODE_TEST compares the orders of accuracy of several explicit Runge-Kutta
%   methods. The non-adaptive ODE solvers are tested on a problem used in
%   E. Hairer, S.P. Norsett, and G. Wanner, Solving Ordinary Differential
%   Equations I, Nonstiff Problems, Springer-Verlag, 1987. 
%   The errors of numerical solutions obtained with several step sizes is
%   plotted against the step size used. For a given solver, the slope of that 
%   line corresponds to the order of the integration method used.  
%

solvers = {'ode1','ode2','ode3','ode4','ode5'};
Nsteps = [800 400 200 100 75 50 20 11];

x0 = 0;
y0 = [1;1;1;1];
xend = 1;
h = (xend-x0)./Nsteps;

err = zeros(length(solvers),length(h));

for i = 1:length(solvers)  
  solver = solvers{i};
  for j = 1:length(Nsteps)
    N = Nsteps(j);
    x = linspace(x0,xend,N+1);
    y = feval(solver,@f,x,y0);
    err(i,j) = max(max(abs(y - yexact(x))));
  end
end  

figure
loglog(h,err);
title('Different slopes for methods of different orders...');
xlabel('log(h)');
ylabel('log(err)');
legend(solvers{:},4);

% -------------------------
function yp = f(x,y)
yp = [  2*x*y(4)*y(1)
       10*x*y(4)*y(1)^5
        2*x*y(4)
       -2*x*(y(3)-1) ];

% -------------------------
function y = yexact(x)
x2 = x(:).^2;
y = zeros(length(x),4);
y(:,1) = exp(sin(x2));
y(:,2) = exp(5*sin(x2));
y(:,3) = sin(x2)+1;
y(:,4) = cos(x2);

:open:

function ode_check_error
%ODE_TEST  Run non-adaptive ODE solvers of different orders.
%   ODE_TEST compares the orders of accuracy of several explicit Runge-Kutta
%   methods. The non-adaptive ODE solvers are tested on a problem used in
%   E. Hairer, S.P. Norsett, and G. Wanner, Solving Ordinary Differential
%   Equations I, Nonstiff Problems, Springer-Verlag, 1987. 
%   The errors of numerical solutions obtained with several step sizes is
%   plotted against the step size used. For a given solver, the slope of that 
%   line corresponds to the order of the integration method used.  
%

solvers = { 'ode1', 'ode2', 'ode3', 'ode4', 'ode5' };
Nsteps  = [ 10 20 40 80 160 320 640 1280 ];

x0   = 0 ;
y0   = [1;1;1;1] ;
xend = 1 ;
h    = (xend-x0)./Nsteps ;

err = zeros(length(solvers),length(h));

for i = 1:length(solvers)
  solver = solvers{i};
  for j = 1:length(Nsteps)
    N = Nsteps(j);
    x = linspace(x0,xend,N+1);
    y = feval(solver,@f,x,y0);
    err(i,j) = max(max(abs(y - yexact(x))));
  end
end  

for i = 1:length(solvers)
  disp( ' ' )
  disp( ['estimate order for scheme = ' solvers{i}] );
  disp( sprintf( 'Error = %-15g', err(i,1) ) ) ;
  for j = 2:length(Nsteps)
    order = log( err(i,j-1)/err(i,j) ) / log( h(j-1)/h(j) ) ;
    disp( sprintf( 'Error = %-15g Order = %-15g', err(i,j), order) ) ;
  end
end  

% -------------------------
function yp = f(x,y)
yp = [  2*x*y(4)*y(1)
       10*x*y(4)*y(1)^5
        2*x*y(4)
       -2*x*(y(3)-1) ];

% -------------------------
function y = yexact(x)
x2 = x.^2;
y  = zeros(length(x),4);
y(:,1) = exp(sin(x2));
y(:,2) = exp(5*sin(x2));
y(:,3) = sin(x2)+1;
y(:,4) = cos(x2);

MATLAB lessons N.2¶