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 ODEy'(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.