ODE45 Numerical Solution of ODEs

Section: Numerical Methods

Usage

function [t,y] = ode45(f,tspan,y0,options,varargin) function SOL = ode45(f,tspan,y0,options,varargin) ode45 is a solver for ordinary differential equations and initial value problems. To solve the ODE
      y'(t) =  f(t,y)
      y(0)  =  y0

over the interval tspan=[t0 t1], you can use ode45. For example, to solve the ode y' = y y(0) = 1 whose exact solution is y(t)=exp(t), over the interval t0=0, t1=3, do

-->       [t,y]=ode45(@(t,y) y,[0 3],1)
Warning: Newly defined variable error shadows a function of the same name.  Use clear error to recover access to the function

k = 
 2 


y = 
    1.0000    1.0030 


t = 

   1.0e-03 * 
         0    3.0000 


k = 
 3 


y = 
    1.0000    1.0030    1.0182 


t = 

   1.0e-02 * 
         0    0.3000    1.8000 


k = 
 4 


y = 
    1.0000    1.0030    1.0182    1.0975 


t = 

   1.0e-02 * 
         0    0.3000    1.8000    9.3000 


k = 
 5 


y = 
    1.0000    1.0030    1.0182    1.0975    1.4814 


t = 
         0    0.0030    0.0180    0.0930    0.3930 


k = 
 6 


y = 
    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997 


t = 
         0    0.0030    0.0180    0.0930    0.3930    0.6930 


k = 
 7 


y = 
    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993 


t = 
         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930 


k = 
 8 


y = 
    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 


t = 
         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 


k = 
 9 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 9

    4.9185 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 9

    1.5930 


k = 
 10 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 10

    4.9185    6.6392 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 10

    1.5930    1.8930 


k = 
 11 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 11

    4.9185    6.6392    8.9620 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 11

    1.5930    1.8930    2.1930 


k = 
 12 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 12

    4.9185    6.6392    8.9620   12.0975 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 12

    1.5930    1.8930    2.1930    2.4930 


k = 
 13 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 13

    4.9185    6.6392    8.9620   12.0975   16.3299 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 13

    1.5930    1.8930    2.1930    2.4930    2.7930 


k = 
 14 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 14

    4.9185    6.6392    8.9620   12.0975   16.3299   20.0854 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 14

    1.5930    1.8930    2.1930    2.4930    2.7930    3.0000 

t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 14

    1.5930    1.8930    2.1930    2.4930    2.7930    3.0000 

y = 
    1.0000 
    1.0030 
    1.0182 
    1.0975 
    1.4814 
    1.9997 
    2.6993 
    3.6437 
    4.9185 
    6.6392 
    8.9620 
   12.0975 
   16.3299 
   20.0854 

If you want a dense output (i.e., an output that also contains an interpolating spline), use instead

-->       SOL=ode45(@(t,y) y,[0 3],1)
Warning: Newly defined variable error shadows a function of the same name.  Use clear error to recover access to the function

k = 
 2 


y = 
    1.0000    1.0030 


t = 

   1.0e-03 * 
         0    3.0000 


k = 
 3 


y = 
    1.0000    1.0030    1.0182 


t = 

   1.0e-02 * 
         0    0.3000    1.8000 


k = 
 4 


y = 
    1.0000    1.0030    1.0182    1.0975 


t = 

   1.0e-02 * 
         0    0.3000    1.8000    9.3000 


k = 
 5 


y = 
    1.0000    1.0030    1.0182    1.0975    1.4814 


t = 
         0    0.0030    0.0180    0.0930    0.3930 


k = 
 6 


y = 
    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997 


t = 
         0    0.0030    0.0180    0.0930    0.3930    0.6930 


k = 
 7 


y = 
    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993 


t = 
         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930 


k = 
 8 


y = 
    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 


t = 
         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 


k = 
 9 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 9

    4.9185 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 9

    1.5930 


k = 
 10 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 10

    4.9185    6.6392 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 10

    1.5930    1.8930 


k = 
 11 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 11

    4.9185    6.6392    8.9620 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 11

    1.5930    1.8930    2.1930 


k = 
 12 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 12

    4.9185    6.6392    8.9620   12.0975 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 12

    1.5930    1.8930    2.1930    2.4930 


k = 
 13 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 13

    4.9185    6.6392    8.9620   12.0975   16.3299 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 13

    1.5930    1.8930    2.1930    2.4930    2.7930 


k = 
 14 


y = 

 Columns 1 to 8

    1.0000    1.0030    1.0182    1.0975    1.4814    1.9997    2.6993    3.6437 

 Columns 9 to 14

    4.9185    6.6392    8.9620   12.0975   16.3299   20.0854 


t = 

 Columns 1 to 8

         0    0.0030    0.0180    0.0930    0.3930    0.6930    0.9930    1.2930 

 Columns 9 to 14

    1.5930    1.8930    2.1930    2.4930    2.7930    3.0000 


SOL = 
    x: 1 14 double array
    y: 1 14 double array
    xe: 
    ye: 
    ie: 
    solver: generic_ode_solver
    interpolant: 1 1 functionpointer array
    idata: 1 1 struct array

You can view the result using

      plot(0:0.01:3,deval(SOL,0:0.01:3))

You will notice that this function is available for "every" value of t, while plot(t,y,'o-') is only available at a few points. The optional argument 'options' is a structure. It may contain any of the following fields: 'AbsTol' - Absolute tolerance, default is 1e-6. 'RelTol' - Relative tolerance, default is 1e-3. 'MaxStep' - Maximum step size, default is (tspan(2)-tspan(1))/10 'InitialStep' - Initial step size, default is maxstep/100 'Stepper' - To override the default Fehlberg integrator 'Events' - To provide an event function 'Projection' - To provide a projection function The varargin is ignored by this function, but is passed to all your callbacks, i.e., f, the event function and the projection function. ==Event Function== The event function can be used to detect situations where the integrator should stop, possibly because the right-hand-side has changed, because of a collision, etc... An event function should look like function [val,isterminal,direction]=event(t,y,...) The return values are: val - the value of the event function. isterminal - whether or not this event should cause termination of the integrator. direction - 1=upcrossings only matter, -1=downcrossings only, 0=both. == Projection function == For geometric integration, you can provide a projection function which will be called after each time step. The projection function has the following signature: function yn=project(t,yn,...); If the output yn is very different from the input yn, the quality of interpolation may decrease.