General form
[T,Y]=solver('F',tspan,y0,options,p1,p2,...) calls F(t,y,flag,p1,p2,...) An example of a stiff problem ============================= Van der Pol with mu=1000 file: vandp.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dy=vandp(t,y) global mu dy=[y(2);mu*(1-y(1)^2)*y(2)-y(1)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global mu mu=1000 %makes the problem stiff tspan=[0 1000]; init=[2;0]; [Trk,Yrk]=ode45('vandp',tspan,init); plot(t,Y(:,1)) title('Sol. of a stiff equ by a non-stiff solver') xlabel('time t'); ylabel('solution Y(:,1)'); >> [Trk,Yrk]=ode45('vandp',tspan,init); DANGER: Memory threshold has been reached. Future memory allocation may result in termination of this process. ??? Out of memory. Type HELP MEMORY for your options. Well, let's quit and start again: global mu mu=1000 %makes the problem stiff tspan=[0 1000]; init=[2;0]; [T,Y]=ode15s('vandp',tspan,init); plot(T,Y(:,1)) title('Sol. of a stiff equ by a stiff solver (ode15s)') xlabel('time t'); ylabel('solution Y(:,1)');
Obtaining sol at specific time points ====================================== The solvers have their own stepsize strategy. Thus, to obtain solutions at other points, you can perform interpolation (by splines). In Matlab V you don't need to do it, though. Instead you can give the timespan as a vector of points. This doesn't affect the computation (at least much) =================================
More extensive ODE-files
Let's modify our vandp.m file: vandp.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [out1,out2,out3]=vandp1(t,y,flag) mu=1; % make mu local and make the problem non-stiff if nargin < 3 | isempty(flag) % return F(t,y), as usual out1= [y(2);mu*(1-y(1)^2)*y(2)-y(1)]; elsif strcmp(flag,'init') % return [tspan,y0,options] out1=[0;20]; out2=[2;0]; out3=odeset('RelTol',1e-4); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Additional parameters are also possible
function [out1,out2,out3]=vandp2(t,y,flag,mu) %Put mu into the parameter list if nargin < 4 | isempty(mu) mu=1; % make mu=1 the default if nargin < 3 | isempty(flag) % return F(t,y), as usual out1= [y(2);mu*(1-y(1)^2)*y(2)-y(1)]; elsif strcmp(flag,'init') % return [tspan,y0,options] out1=[0;20]; out2=[2;0]; out3=odeset('RelTol',1e-4); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The nr. of additional params is arbitrary.Examples
>> type rigidode function [out1,out2,out3] = rigidode(t,y,flag) %RIGIDODE Euler equations of a rigid body without external forces. % A standard test problem for non-stiff solvers proposed by Krogh. The % analytical solutions are Jacobian elliptic functions accessible in % MATLAB. The interval here is about 1.5 periods; it is that for which % solutions are plotted on p. 243 of Shampine and Gordon. % % RIGIDODE([],[],'init') returns the default TSPAN, Y0, and OPTIONS values % for this problem. These values are retrieved by an ODE Suite solver if % the solver is invoked with empty TSPAN or Y0 arguments. This example % does not set any OPTIONS, so the third output argument is set to empty % [] instead of an OPTIONS structure created with ODESET. % % L. F. Shampine and M. K. Gordon, Computer Solution of Ordinary % Differential Equations, W.H. Freeman & Co., 1975. % % See also ODE45, ODE23, ODE113, ODESET, ODEFILE. % Mark W. Reichelt and Lawrence F. Shampine, 3-23-94, 4-19-94 % Copyright (c) 1984-96 by The MathWorks, Inc. % $Revision: 1.5 $ $Date: 1996/06/26 16:23:46 $ if nargin < 3 | isempty(flag) % Return dy/dt = f(t,y). out1 = [y(2)*y(3); -y(1)*y(3); -0.51*y(1)*y(2)]; else switch(flag) case 'init' % Used only if TSPAN or Y0 is empty. % Return default TSPAN, Y0, and OPTIONS. out1 = [0; 12]; out2 = [0; 1; 1]; out3 = []; otherwise error(['Unknown flag ''' flag '''.']); end endExample calls
1. rigidode
>> rigidode([],[],'init') ans = 0 12 >> [o1,o2]=rigidode([],[],'init') o1 = 0 12 o2 = 0 1 1 >> ode45('rigidode') % Call with no arguments % Remember, the odeplot-function is the default output function.
vdpode
function [out1,out2,out3] = vdpode(t,y,flag,mu) %VDPODE Parameterizable van der Pol equation (stiff for large mu). % VDPODE(T,Y) or VDPODE(T,Y,[],MU) returns the derivatives vector for the % van der Pol equation. By default, MU is 1, and the problem is not % stiff. Optionally, pass in the MU parameter as an additional parameter % to an ODE Suite solver. The problem becomes more stiff as MU is % increased. % % When MU is 1000 the equation is in relaxation oscillation, and the % problem becomes very stiff. The limit cycle has portions where the % solution components change slowly and the problem is quite stiff, % alternating with regions of very sharp change where it is not stiff % (quasi-discontinuities). The initial conditions are close to an area of % slow change so as to test schemes for the selection of the initial step % size. % % VDPODE(T,Y,'jacobian') or VDPODE(T,Y,'jacobian',MU) returns the Jacobian % matrix dF/dY evaluated analytically at (T,Y). By default, the stiff % solvers of the ODE Suite approximate Jacobian matrices numerically. % However, if the ODE Solver property Jacobian is set to 'on' with ODESET, % a solver calls the ODE file with the flag 'jacobian' to obtain dF/dY. % Providing the solvers with an analytic Jacobian is not necessary, but it % can improve the reliability and efficiency of integration. % % VDPODE([],[],'init') returns the default TSPAN, Y0, and OPTIONS values % for this problem (see RIGIDODE). The ODE solver property Vectorized is % set to 'on' with ODESET because VDPODE is coded so that calling % VDPODE(T,[Y1 Y2 ...]) returns [VDPODE(T,Y1) VDPODE(T,Y2) ...] for % scalar time T and vectors Y1,Y2,... The stiff solvers of the ODE Suite % take advantage of this feature when approximating the columns of the % Jacobian numerically. % % L. F. Shampine, Evaluation of a test set for stiff ODE solvers, ACM % Trans. Math. Soft., 7 (1981) pp. 409-420. % % See also ODE15S, ODE23S, ODESET, ODEFILE. % Mark W. Reichelt and Lawrence F. Shampine, 3-23-94, 4-19-94 % Copyright (c) 1984-96 by The MathWorks, Inc. % $Revision: 1.7 $ $Date: 1996/07/08 12:31:28 $ if nargin < 4 | isempty(mu) mu = 1; end if nargin < 3 | isempty(flag) % Return dy/dt = f(t,y). out1 = [y(2,:); (mu*(1-y(1,:).^2).*y(2,:) - y(1,:))]; % vectorized else switch(flag) case 'init' % Used only if TSPAN or Y0 is empty. % Return default TSPAN, Y0, and OPTIONS. out1 = [0; max(20,3*mu)]; % several periods out2 = [2; 0]; out3 = odeset('Vectorized','on'); case 'jacobian' % Used only if odeset('Jacobian','on'). % Return Jacobian matrix J(t,y) = df/dy, evaluating it analytically. out1 = [ 0 1 (-2*mu*y(1)*y(2) - 1) (mu*(1-y(1)^2)) ]; otherwise error(['Unknown flag ''' flag '''.']); end endballode
function [out1,out2,out3] = ballode(t,y,flag) %BALLODE Equations of motion for a bouncing ball, used by BALLDEM. % BALLODE(T,Y) returns the derivatives vector for the equations of motion % for a bouncing ball. This ODE file can be used to test the % zero-crossing location capabilities of the ODE Suite solvers. % % BALLODE(T,Y,'events') returns a zero-crossing vector VALUE evaluated at % (T,Y) and two constant vectors ISTERMINAL and DIRECTION. By default, % the solvers of the ODE Suite do not locate zero-crossings. However, if % the ODE solver property Events is set to 'on' with ODESET, the solver % calls the ODE file with the flag 'events'. The ODE file returns the % information that the solver uses to locate zero-crossings of the % elements in the VALUE vector. The VALUE vector may be of any length. % It is evaluated at the beginning and end of a step, and if any elements % make transitions to, from, or through zero (with the directionality % specified in DIRECTION), then the zero-crossing point is located. The % ISTERMINAL vector consists of logical 1's and 0's, enabling you to % specify whether or not a zero-crossing of the corresponding VALUE % element halts the integration. The DIRECTION vector enables you to % specify a desired directionality, positive (1), negative (-1), and don't % care (0) for each VALUE element. % % BALLODE also responds to the flag 'init' (see RIGIDODE). % % See also BALLDEM, ODE45, ODESET, ODEFILE. % Mark W. Reichelt and Lawrence F. Shampine, 1/3/95 % Copyright (c) 1984-96 by The MathWorks, Inc. % $Revision: 1.7 $ $Date: 1996/07/08 12:29:29 $ if nargin < 3 | isempty(flag) % Return dy/dt = f(t,y). out1 = [y(2); -9.8]; else switch(flag) case 'init' % Used only if TSPAN or Y0 is empty. % Return default TSPAN, Y0, and OPTIONS. out1 = [0; 10]; out2 = [0; 20]; out3 = odeset('Events','on'); case 'events' % Used only if odeset('Events','on'). % Return event vectors VALUE, ISTERMINAL, and DIRECTION. % Locate zero-crossings of both height and velocity. out1 = y; % [height; velocity] out2 = [1; 0]; % stop at zeros of height % event only when value(1) is decreasing, don't care for value(2) out3 = [-1; 0]; otherwise error(['Unknown flag ''' flag '''.']); end end Driver function for ballode >> type balldem function [ttout,yout,teout,yeout,ieout] = balldem %BALLDEM Driver function for the ODE Suite bouncing ball simulation. % BALLDEM is an example of repeated event location, where the initial % conditions are changed after each terminal event. The BALLDEM function % computes ten bounces of the bouncing ball problem, BALLODE with event % functions. Each bounce is solved with a call to ODE45. The speed of % the ball is attenuated by 0.9 after each bounce, and the simulation is % stopped when the time for a bounce has decreased to a minimum length. % Note that the event function also locates the peak of each bounce. % % [T,Y,TE,YE,IE] = BALLDEM returns the accumulation of the trajectories % (T,Y) and the event location times TE, solution values YE, and event % types IE. % % See also BALLODE, ODE45. % Mark W. Reichelt and Lawrence F. Shampine, 1/3/95 % Copyright (c) 1984-96 by The MathWorks, Inc. % $Revision: 1.7 $ $Date: 1996/10/02 14:00:34 $ tstart = 0; tfinal = 30; y0 = [0; 20]; refine = 4; options = odeset('Events','on','OutputFcn','odeplot','OutputSel',1,... 'Refine',refine); clf reset % deletes any stop button set(gca,'xlim',[0 30],'ylim',[0 25]); box on hold on; tout = tstart; yout = y0.'; teout = []; yeout = []; ieout = []; for i = 1:10 % Solve until the first terminal event. [t,y,te,ye,ie] = ode23('ballode',[tstart tfinal],y0,options); if ~ishold hold on end % Accumulate output. nt = length(t); tout = [tout; t(2:nt)]; yout = [yout; y(2:nt,:)]; teout = [teout; te]; % Events at tstart are never reported. yeout = [yeout; ye]; ieout = [ieout; ie]; ud = get(gcf,'UserData'); if ud.stop break; end % Set the new initial conditions, with .9 attenuation. y0(1) = 0; y0(2) = -.9*y(nt,2); % A good guess of a valid first timestep is the length of the last valid % timestep, so use it for faster computation. 'refine' is 4 by default. options = odeset(options,'InitialStep',t(nt)-t(nt-refine),... 'MaxStep',t(nt)-t(1)); tstart = t(nt); end hold off if nargout > 0 ttout = tout; end >> ballode