matlab中ode45函数编写.docx

1、matlab中ode45函数编写function varargout = ode45(ode,tspan,y0,options,varargin)%ODE45 Solve non-stiff differential equations, medium order method.% TOUT,YOUT = ODE45(ODEFUN,TSPAN,Y0) with TSPAN = T0 TFINAL integrates % the system of differential equations y = f(t,y) from time T0 to TFINAL % with initial c

2、onditions Y0. ODEFUN is a function handle. For a scalar T% and a vector Y, ODEFUN(T,Y) must return a column vector corresponding % to f(t,y). Each row in the solution array YOUT corresponds to a time % returned in the column vector TOUT. To obtain solutions at specific % times T0,T1,.,TFINAL (all in

3、creasing or all decreasing), use TSPAN = % T0 T1 . TFINAL. % % TOUT,YOUT = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default% integration properties replaced by values in OPTIONS, an argument created% with the ODESET function. See ODESET for details. Commonly used options % are scalar rela

4、tive error tolerance RelTol (1e-3 by default) and vector% of absolute error tolerances AbsTol (all components 1e-6 by default).% If certain components of the solution must be non-negative, use% ODESET to set the NonNegative property to the indices of these% components.% % ODE45 can solve problems M(

5、t,y)*y = f(t,y) with mass matrix M that is% nonsingular. Use ODESET to set the Mass property to a function handle % MASS if MASS(T,Y) returns the value of the mass matrix. If the mass matrix % is constant, the matrix can be used as the value of the Mass option. If% the mass matrix does not depend on

6、 the state variable Y and the function% MASS is to be called with one input argument T, set MStateDependence to% none. ODE15S and ODE23T can solve problems with singular mass matrices. % TOUT,YOUT,TE,YE,IE = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) with the Events% property in OPTIONS set to a function handle

7、 EVENTS, solves as above % while also finding where functions of (T,Y), called event functions, % are zero. For each function you specify whether the integration is % to terminate at a zero and whether the direction of the zero crossing % matters. These are the three column vectors returned by EVENT

8、S: % VALUE,ISTERMINAL,DIRECTION = EVENTS(T,Y). For the I-th event function: % VALUE(I) is the value of the function, ISTERMINAL(I)=1 if the integration % is to terminate at a zero of this event function and 0 otherwise. % DIRECTION(I)=0 if all zeros are to be computed (the default), +1 if only % zer

9、os where the event function is increasing, and -1 if only zeros where % the event function is decreasing. Output TE is a column vector of times % at which events occur. Rows of YE are the corresponding solutions, and % indices in vector IE specify which event occurred. % SOL = ODE45(ODEFUN,T0 TFINAL

10、,Y0.) returns a structure that can be% used with DEVAL to evaluate the solution or its first derivative at % any point between T0 and TFINAL. The steps chosen by ODE45 are returned % in a row vector SOL.x. For each I, the column SOL.y(:,I) contains % the solution at SOL.x(I). If events were detected

11、, SOL.xe is a row vector % of points at which events occurred. Columns of SOL.ye are the corresponding % solutions, and indices in vector SOL.ie specify which event occurred. % Example % t,y=ode45(vdp1,0 20,2 0); % plot(t,y(:,1);% solves the system y = vdp1(t,y), using the default relative error% to

12、lerance 1e-3 and the default absolute tolerance of 1e-6 for each% component, and plots the first component of the solution. % % Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y):% float: double, single% See also % other ODE solvers: ODE23, ODE113, ODE15S, ODE23S, ODE23T, ODE23TB % im

13、plicit ODEs: ODE15I% options handling: ODESET, ODEGET% output functions: ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT% evaluating solution: DEVAL% ODE examples: RIGIDODE, BALLODE, ORBITODE% function handles: FUNCTION_HANDLE % NOTE: % The interpretation of the first input argument of the ODE solvers and % s

14、ome properties available through ODESET have changed in MATLAB 6.0.% Although we still support the v5 syntax, any new functionality is % available only with the new syntax. To see the v5 help, type in % the command line % more on, type ode45, more off % NOTE:% This portion describes the v5 syntax of

15、 ODE45.% T,Y = ODE45(F,TSPAN,Y0) with TSPAN = T0 TFINAL integrates the% system of differential equations y = F(t,y) from time T0 to TFINAL with% initial conditions Y0. F is a string containing the name of an ODE% file. Function F(T,Y) must return a column vector. Each row in% solution array Y corres

16、ponds to a time returned in column vector T. To% obtain solutions at specific times T0, T1, ., TFINAL (all increasing% or all decreasing), use TSPAN = T0 T1 . TFINAL.% % T,Y = ODE45(F,TSPAN,Y0,OPTIONS) solves as above with default% integration parameters replaced by values in OPTIONS, an argument% c

17、reated with the ODESET function. See ODESET for details. Commonly% used options are scalar relative error tolerance RelTol (1e-3 by% default) and vector of absolute error tolerances AbsTol (all% components 1e-6 by default).% % T,Y = ODE45(F,TSPAN,Y0,OPTIONS,P1,P2,.) passes the additional% parameters

18、 P1,P2,. to the ODE file as F(T,Y,FLAG,P1,P2,.) (see% ODEFILE). Use OPTIONS = as a place holder if no options are set.% % It is possible to specify TSPAN, Y0 and OPTIONS in the ODE file (see% ODEFILE). If TSPAN or Y0 is empty, then ODE45 calls the ODE file% TSPAN,Y0,OPTIONS = F(,init) to obtain any

19、values not supplied% in the ODE45 argument list. Empty arguments at the end of the call list% may be omitted, e.g. ODE45(F).% % ODE45 can solve problems M(t,y)*y = F(t,y) with a mass matrix M that is% nonsingular. Use ODESET to set Mass to M, M(t), or M(t,y) if the% ODE file is coded so that F(T,Y,m

20、ass) returns a constant,% time-dependent, or time- and state-dependent mass matrix, respectively.% The default value of Mass is none. ODE15S and ODE23T can solve problems% with singular mass matrices. % % T,Y,TE,YE,IE = ODE45(F,TSPAN,Y0,OPTIONS) with the Events property in% OPTIONS set to on, solves

21、 as above while also locating zero crossings% of an event function defined in the ODE file. The ODE file must be% coded so that F(T,Y,events) returns appropriate information. See% ODEFILE for details. Output TE is a column vector of times at which% events occur, rows of YE are the corresponding solu

22、tions, and indices in% vector IE specify which event occurred.% % See also ODEFILE % ODE45 is an implementation of the explicit Runge-Kutta (4,5) pair of% Dormand and Prince called variously RK5(4)7FM, DOPRI5, DP(4,5) and DP54.% It uses a free interpolant of order 4 communicated privately by% Dorman

23、d and Prince. Local extrapolation is done. % Details are to be found in The MATLAB ODE Suite, L. F. Shampine and% M. W. Reichelt, SIAM Journal on Scientific Computing, 18-1, 1997. % Mark W. Reichelt and Lawrence F. Shampine, 6-14-94% Copyright 1984-2009 The MathWorks, Inc.% $Revision: 5.74.4.10 $ $D

24、ate: 2009/04/21 03:24:15 $ solver_name = ode45; % Check inputsif nargin 4 options = ; if nargin 3 y0 = ; if nargin 2 tspan = ; if nargin 0); % t,y,. = odeXX(.)% There might be no output requested. sol = ; f3d = ; if output_sol sol.solver = solver_name; sol.extdata.odefun = ode; sol.extdata.options =

25、 options; sol.extdata.varargin = varargin; end % Handle solver argumentsneq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, . options, threshold, rtol, normcontrol, normy, hmax, htry, htspan, dataType = . odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);n

26、fevals = nfevals + 1; % Handle the outputif nargout 0 outputFcn = odeget(options,OutputFcn,fast);else outputFcn = odeget(options,OutputFcn,odeplot,fast);endoutputArgs = ; if isempty(outputFcn) haveOutputFcn = false;else haveOutputFcn = true; outputs = odeget(options,OutputSel,1:neq,fast); if isa(out

27、putFcn,function_handle) % With MATLAB 6 syntax pass additional input arguments to outputFcn. outputArgs = varargin; end endrefine = max(1,odeget(options,Refine,4,fast);if ntspan 2 outputAt = RequestedPoints; % output only at tspan pointselseif refine 0 % non-trivial mass matrix Msingular = odeget(op

28、tions,MassSingular,no,fast); if strcmp(Msingular,maybe) warning(MATLAB:ode45:MassSingularAssumedNo,ODE45 assumes . MassSingular is no. See ODE15S or ODE23T.); elseif strcmp(Msingular,yes) error(MATLAB:ode45:MassSingularYes,. MassSingular cannot be yes for this solver. See ODE15S . or ODE23T.); end % Incorporate the mass matrix into odeFcn and odeArgs. odeFcn,odeArgs = odemassexplicit(FcnHandlesUsed,Mtype,odeFcn,odeArgs,Mfun,M); f0