Matlab Demo

More on Inline Functions

As I was preparing these notes, I learned something new about inline functions.

A standard inline function stores away the definition (as a string), and does not evaluate it until the function is called. At that point, any variable names other than the arguments are undefined: 

>> a = 2
a =
     2
>> f = inline('a*x','x')
f =
     Inline function:
     f(x) = a*x
>> f(5)
??? Error using ==> inlineeval
Undefined function or variable 'a'.

An inline function defined with the other notation evaluates variables that are not arguments. I did not see that mentioned in the manual: 

>> a = 2
a =
     2
>> f = @(x) a*x
f = 
    @(x) a*x
>> f(5)
ans =
    10
>> a = 3
a =
     3
>> f(5)
ans =
    10

Numerical Differentiation

First Derivative

The only differentiation routines built into Matlab are ppder (for polynomials), fnder (for splines), and diff (symbolic differentiation). For general functions you need to program your own formulas.

Example: Estimate the derivative of ex at x=1.

The standard forward / backward / centered difference formulas produce

>> h = 0.1; x = 1;
>> forward = (exp(x+h) - exp(x))/h, error = forward - exp(1)
forward =
   2.858841954873879
error =
   0.140560126414833
>> backward = (exp(x) - exp(x - h))/h, error = backward - exp(1)
backward =
   2.586787173020957
error =
  -0.131494655438089
>> centered = (exp(x+h) - exp(x-h))/(2*h), error = centered - exp(1)
centered =
   2.722814563947418
error =
   0.0045327354883

If you have irregular given data points, you can interpolate and differentiate the interpolant, using either polynomials or splines. 

>> xx = [0.5,1.2,2,3]; yy = exp(xx);
>> p = polyfit(xx,yy,3);
>> p_derivative = polyval(polyder(p),x)
p_derivative =
   2.538765895391738
>> error = p_derivative - exp(1)
error =
  -0.179515933067308
>> s = spline(xx,yy);
>> s_derivative = ppval(fnder(s),x)
s_derivative =
   2.538765895391735
>> error = s_derivative - exp(1)
error =
  -0.179515933067310

Second Derivative

>> centered = (exp(x+h) - 2*exp(x) + exp(x-h))/h^2
centered =
   2.720547818529217
>> error = centered - exp(1)
error =
   0.002265990070172
>> s_derivative = ppval(fnder(s,2),x)
s_derivative =
   2.237730774428740
>> error = s_derivative - exp(1)
error =
  -0.480551054030306
>> p_derivative = polyval(polyder(polyder(p)),x)
p_derivative =
   2.237730774428767
>> error = p_derivative - exp(1)
error =
  -0.480551054030279

Numerical Integration

One-Dimensional Integral

The standard integration routines in Matlab are quad and quadl.

Example: Integrate ex from 0 to 1. 

>> [I,neval] = quad('exp',0,1), error = I - (exp(1)-1)
I =
   1.718281832201628
neval =
    13
error =
    3.742582688204266e-009
>> [I,neval] = quad('exp',0,1,1.e-12), error = I - (exp(1)-1)
I =
   1.718281828459045
neval =
   189
error =
   -2.220446049250313e-016
>> [I,neval] = quadl('exp',0,1), error = I - (exp(1)-1)
I =
   1.718281828459183
neval =
    18
error =
    1.374456104485944e-013
>> [I,neval] = quadl('exp',0,1,1.e-12), error = I - (exp(1)-1)
I =
   1.718281828459183
neval =
    18
error =
    1.374456104485944e-013

In my opinion, numerical integration is a weak spot in Matlab. There ought to be some routines for doing mildly singular integrals, such as sqrt(x), or more seriously singular integrals, such as 1/sqrt(x). Notice that integrating sqrt(x) takes a lot more evaluations than the exponential function. Most of them are near 0, where the singularity is. 

>> [I,neval] = quadl('sqrt',0,1), error = I - 2/3
I =
   0.666665345816982
neval =
    78
error =
   -1.320849684627312e-006
>> [I,neval] = quadl('sqrt',0,1,1.e-12), error = I - 2/3
I =
   0.666666666665878
neval =
   588
error =
   -7.881473251813986e-013

Multidimensional Integrals

Example: Integrate ex+y sin x over the square [0,1] x [1,2]. Let's do the exact integral first:

>> true = int(int(sym('exp(x+y)*sin(x)'),'x',0,1),'y',1,2)
true =
-1/2*exp(3)*cos(1)+1/2*exp(3)*sin(1)+1/2*exp(2)+1/2*exp(2)*cos(1)-1/2*exp(2)*sin(1)-1/2*exp(1)
>> true = double(true)
true =
   4.247278313748523

Method 1: Do an interated integral.

This is considerably more tricky than it seems. I tried to do it with inline functions and failed. You have to use functions stored in files. Defining the integrand is easy: 

function z = f(x,y)
z = exp(x + y) .* sin(x);

The inner integral is

function z = inner(y)
% do x integration for fixed y
z = zeros(size(y));
for i = 1:length(z)
    z(i) = quad('f',0,1,[],[],y(i));
end

Notice some things here:

  • We are using an undocumented feature of quad, where any arguments after the first 5 are passed on the function f. That lets us get the y argument in there.
  • The code 
    z = quad('f',0,1,[],[],y)
    does not work. We have to put in an explicit loop and pass each y(i) one by one.
  • It has to be x integration first, then y.

After the setup it works quite well: 

>> I = quad('inner',1,2), error = I - true
I =
   4.247278308991272
error =
   -4.757250593456774e-009

Method 2: Use dblquad:

>> I = dblquad('f',0,1,1,2), error = I - true
I =
   4.247278308991272
error =
   -4.757250593456774e-009

The fact that you get exactly the same number is no coincidence: behind the scenes, dblquad does precisely what we did in Method 1.

The only problem with dblquad is that it requires a rectangular area. If you want to integrate over a different area, you need to write your function so that it has value 0 outside the area, and the resulting discontinuity will slow it down a bit.

Do-It-Yourself ODE Solvers

Example: y'(t) = -2 t y2(t), y(0) = 1, true solution y(t) = 1/(1 + t2). Use step size 0.1 from 0 to 2.

Euler's Method:

>> f = inline('-2*t.*y.^2','t','y'); y0 = 1;
>> true = inline('-1./(1 + t.^2)');
>> h = 0.1; t = 0:h:2; y_true = true(t);
>> y_euler = zeros(size(t)); y_euler(1) = y0;
>> for i = 1:length(t)-1
       y_euler(i+1) = y_euler(i) + h * f(t(i),y_euler(i));
   end
>> plot(t,y_true,'k',t,y_euler,'r')
>> title('Euler''s method')

Euler's method

Runge-Kutta:

>> y_rk = zeros(size(t)); y_rk(1) = y0;
>> for i = 1:length(t)-1
       k1 = h*f(t(i),y_rk(i));
       k2 = h*f(t(i)+0.5*h,y_rk(i)+0.5*k1);
       k3 = h*f(t(i)+0.5*h,y_rk(i)+0.5*k2);
       k4 = h*f(t(i)+h,y_rk(i)+k3);
       y_rk(i+1) = y_rk(i) + (k1 + 2*k2 + 2*k3 + k4)/6;
   end
>> plot(t,y_true,'k',t,y_rk,'r')
>> title('Runge-Kutta method')

Runge-Kutta

Multistep:

Notice that we need Runge-Kutta to provide start-up values. We also need a separate vector for the f values. 

>> y_abm = zeros(size(t)); y_abm(1:4) = y_rk(1:4);
>> f_abm = zeros(size(t)); f_abm(1:4) = f(t(1:4),y_abm(1:4));
>> for i = 4:length(t)-1
       y_abm(i+1) = y_abm(i) + h/24 * (55*f_abm(i) - 59*f_abm(i-1) + 37*f_abm(i-2) - 9*f_abm(i-3));  % P
       f_abm(i+1) = f(t(i+1),y_abm(i+1))                                                             % E
       y_abm(i+1) = y_abm(i) + h/24 * (9*f_abm(i+1) + 19*f_abm(i) - 5*f_abm(i-1) +f_abm(i-2));       % C
       f_abm(i+1) = f(t(i+1),y_abm(i+1));                                                            % E
   end
>> plot(t,y_true,'k',t,y_abm,'r')
>> title('Adams-Bashforth-Moulton method')

Adams-Bashforth-Moulton

Matlab ODE Solvers

ODE45

We do the same example as before:

>> [tt,y_ode45] = ode45(f,[0,2],y0);
>> plot(t,y_true,'k',tt,y_ode45,'r')
>> title('ODE45')
>> tt
tt =
                   0
   0.050000000000000
   0.100000000000000
   0.150000000000000
   ...
   2.000000000000000
>> s = spline(tt,y_ode45);
>> y2_ode45 = ppval(s,t);

ODE45

Comments:

  • ODE45 chooses its own steps. To get values at the points you want, you have to do spline interpolation afterwards.
  • In this simple example, ODE45 took steps of constant size, but usually the step size is more irregular.

You can set various options with the odeset command. For example, let's do the same calculation with reduced accuracy:

>> odeset
          AbsTol: [ positive scalar or vector {1e-6} ]
          RelTol: [ positive scalar {1e-3} ]
     NormControl: [ on | {off} ]
     NonNegative: [ vector of integers ]
       OutputFcn: [ function_handle ]
       OutputSel: [ vector of integers ]
          Refine: [ positive integer ]
           Stats: [ on | {off} ]
     InitialStep: [ positive scalar ]
         MaxStep: [ positive scalar ]
             BDF: [ on | {off} ]
        MaxOrder: [ 1 | 2 | 3 | 4 | {5} ]
        Jacobian: [ matrix | function_handle ]
        JPattern: [ sparse matrix ]
      Vectorized: [ on | {off} ]
            Mass: [ matrix | function_handle ]
MStateDependence: [ none | {weak} | strong ]
       MvPattern: [ sparse matrix ]
    MassSingular: [ yes | no | {maybe} ]
    InitialSlope: [ vector ]
          Events: [ function_handle ]

>> options = odeset('AbsTol',1.e-3);
>> [tt,y_ode45] = ode45(f,[0,2],y0,options);

In this particular example, that has no effect. You still get the same values for tt.

Higher Order ODEs

As an example of a higher order ODE, or equivalently a system of ODEs, consider the van der Pol equation

y ''(t) - μ [1 - y2(t)] y '(t) + y(t) = 0

Solutions of this equation become approximately periodic with a period of order μ. The time scales involved are of order μ and 1/μ, so for large values of μ this equation becomes stiff.

Runge-Kutta:

All the standard methods are unchanged when you go from scalar equations to vector equations. I will just illustrate that with the Runge-Kutta method here. We use the van der Pol equation with μ=2 (non-stiff). 

>> mu = 2; z0 = [1;0]; h = 0.1;
>> f = @(t,z) [z(2);mu*(1-z(1)^2)*z(2)-z(1)];
>> t_rk = 0:h:20; z_rk = zeros(2,length(t_rk)); z_rk(:,1) = z0;
>> for i = 1:length(t_rk)-1
       k1 = h*f(t_rk(i),z_rk(:,i));
       k2 = h*f(t_rk(i)+0.5*h,z_rk(:,i)+0.5*k1);
       k3 = h*f(t_rk(i)+0.5*h,z_rk(:,i)+0.5*k2);
       k4 = h*f(t_rk(i)+h,z_rk(:,i)+k3);
       z_rk(:,i+1) = z_rk(:,i) + (k1 + 2*k2 + 2*k3 + k4)/6;
   end
>> plot(t_rk,z_rk(1,:))
>> title('\mu=2, Runge-Kutta')

Runge-Kutta

Notice two things here:

  • In the definition of f, we have to use the @(t,z) notation. The inline command does not work, because the value of μ would not be available. Using a function file would also work, of course, but we are trying to avoid that.
  • Instead of z(i), we use z(:,i). We need to pick a column from a matrix, not just an entry from a vector.

ODE45

Passing the value of μ is a bit tricky. We have two choices:

  • Make f into a function file, and pass μ as a global variable.
  • Use the undocumented feature of ODE45 that anything after the standard parameters gets passed along to f as extra parameters.

We use the second method. 

>> f = @(t,z,mu) [z(2);mu*(1-z(1)^2)*z(2)-z(1)];
>> [t,z_ode45] = ode45(f,[0,20],z0,[],mu);
>> size(t)
ans =
   309     1
>> plot(t_rk,z_rk(1,:),'k',t,z_ode45(:,1),'r')
>> title('\mu=2, RK (black) versus ODE45 (red)')

RK versus ODE45

I ran into one problem here: notice that it says z_rk(1,:), but z_ode45(:,1). Even though z0 is a column vector, and f is written to produce column vectors, ODE45 outputs row vectors. Live with it.

Stiff ODEs

This time we use μ=100, which makes the problem quite stiff. You can still solve it with ODE45, but it is much faster to use ODE23s: 

>> mu = 100; z0 = [1;0];
>> f = @(t,z,mu) [z(2);mu*(1-z(1)^2)*z(2)-z(1)];
>> tic;[t_ode45,z_ode45] = ode45(f,[0,400],z0,[],mu);toc;
Elapsed time is 11.050036 seconds.
>> size(t_ode45)
ans =
       87841           1
>> tic;[t_ode23s,z_ode23s] = ode23s(f,[0,400],z0,[],mu);toc;
Elapsed time is 0.918247 seconds.
>> size(t_ode23s)
ans =
   791     1
>> plot(t_ode23s,z_ode23s(:,1),'k',t_ode45,z_ode45(:,1),'r')
>> title('\mu=100, ODE23s (black) versus ODE45 (red)')

ODE23s

For comparison, I tried to run standard Runge-Kutta with step size h=0.01 on this problem. It did OK for a little bit, then this happened:

>> z_rk(1,50)
ans =
  -2.495706117699899
>> z_rk(1,51)
ans =
    3.782800431352553e+005
>> z_rk(1,52)
ans =
    7.518454240602221e+095

 

Last Updated: March 12, 2009