matlab中ode45函数编写.docx

matlab中ode45函数编写.docx

《matlab中ode45函数编写.docx》由会员分享，可在线阅读，更多相关《matlab中ode45函数编写.docx（23页珍藏版）》请在冰豆网上搜索。

matlab中ode45函数编写

functionvarargout=ode45（ode,tspan,y0,options,varargin）

%ODE45Solvenon-stiffdifferentialequations,mediumordermethod.

%[TOUT,YOUT]=ODE45（ODEFUN,TSPAN,Y0）withTSPAN=[T0TFINAL]integrates

%thesystemofdifferentialequationsy'=f（t,y）fromtimeT0toTFINAL

%withinitialconditionsY0.ODEFUNisafunctionhandle.ForascalarT

%andavectorY,ODEFUN（T,Y）mustreturnacolumnvectorcorresponding

%tof（t,y）.EachrowinthesolutionarrayYOUTcorrespondstoatime

%returnedinthecolumnvectorTOUT.Toobtainsolutionsatspecific

%timesT0,T1,...,TFINAL（allincreasingoralldecreasing）,useTSPAN=

%[T0T1...TFINAL].

%

%[TOUT,YOUT]=ODE45（ODEFUN,TSPAN,Y0,OPTIONS）solvesasabovewithdefault

%integrationpropertiesreplacedbyvaluesinOPTIONS,anargumentcreated

%withtheODESETfunction.SeeODESETfordetails.Commonlyusedoptions

%arescalarrelativeerrortolerance'RelTol'（1e-3bydefault）andvector

%ofabsoluteerrortolerances'AbsTol'（allcomponents1e-6bydefault）.

%Ifcertaincomponentsofthesolutionmustbenon-negative,use

%ODESETtosetthe'NonNegative'propertytotheindicesofthese

%components.

%

%ODE45cansolveproblemsM（t,y）*y'=f（t,y）withmassmatrixMthatis

%nonsingular.UseODESETtosetthe'Mass'propertytoafunctionhandle

%MASSifMASS（T,Y）returnsthevalueofthemassmatrix.Ifthemassmatrix

%isconstant,thematrixcanbeusedasthevalueofthe'Mass'option.If

%themassmatrixdoesnotdependonthestatevariableYandthefunction

%MASSistobecalledwithoneinputargumentT,set'MStateDependence'to

%'none'.ODE15SandODE23Tcansolveproblemswithsingularmassmatrices.

%

%[TOUT,YOUT,TE,YE,IE]=ODE45（ODEFUN,TSPAN,Y0,OPTIONS）withthe'Events'

%propertyinOPTIONSsettoafunctionhandleEVENTS,solvesasabove

%whilealsofindingwherefunctionsof（T,Y）,calledeventfunctions,

%arezero.Foreachfunctionyouspecifywhethertheintegrationis

%toterminateatazeroandwhetherthedirectionofthezerocrossing

%matters.ThesearethethreecolumnvectorsreturnedbyEVENTS:

%[VALUE,ISTERMINAL,DIRECTION]=EVENTS（T,Y）.FortheI-theventfunction:

%VALUE（I）isthevalueofthefunction,ISTERMINAL（I）=1iftheintegration

%istoterminateatazeroofthiseventfunctionand0otherwise.

%DIRECTION（I）=0ifallzerosaretobecomputed（thedefault）,+1ifonly

%zeroswheretheeventfunctionisincreasing,and-1ifonlyzeroswhere

%theeventfunctionisdecreasing.OutputTEisacolumnvectoroftimes

%atwhicheventsoccur.RowsofYEarethecorrespondingsolutions,and

%indicesinvectorIEspecifywhicheventoccurred.

%

%SOL=ODE45（ODEFUN,[T0TFINAL],Y0...）returnsastructurethatcanbe

%usedwithDEVALtoevaluatethesolutionoritsfirstderivativeat

%anypointbetweenT0andTFINAL.ThestepschosenbyODE45arereturned

%inarowvectorSOL.x.ForeachI,thecolumnSOL.y（:

I）contains

%thesolutionatSOL.x（I）.Ifeventsweredetected,SOL.xeisarowvector

%ofpointsatwhicheventsoccurred.ColumnsofSOL.yearethecorresponding

%solutions,andindicesinvectorSOL.iespecifywhicheventoccurred.

%

%Example

%[t,y]=ode45（@vdp1,[020],[20]）;

%plot（t,y（:

1））;

%solvesthesystemy'=vdp1（t,y）,usingthedefaultrelativeerror

%tolerance1e-3andthedefaultabsolutetoleranceof1e-6foreach

%component,andplotsthefirstcomponentofthesolution.

%

%ClasssupportforinputsTSPAN,Y0,andtheresultofODEFUN（T,Y）:

%float:

double,single

%

%Seealso

%otherODEsolvers:

ODE23,ODE113,ODE15S,ODE23S,ODE23T,ODE23TB

%implicitODEs:

ODE15I

%optionshandling:

ODESET,ODEGET

%outputfunctions:

ODEPLOT,ODEPHAS2,ODEPHAS3,ODEPRINT

%evaluatingsolution:

DEVAL

%ODEexamples:

RIGIDODE,BALLODE,ORBITODE

%functionhandles:

FUNCTION_HANDLE

%

%NOTE:

%TheinterpretationofthefirstinputargumentoftheODEsolversand

%somepropertiesavailablethroughODESEThavechangedinMATLAB6.0.

%Althoughwestillsupportthev5syntax,anynewfunctionalityis

%availableonlywiththenewsyntax.Toseethev5help,typein

%thecommandline

%moreon,typeode45,moreoff

%NOTE:

%Thisportiondescribesthev5syntaxofODE45.

%

%[T,Y]=ODE45（'F',TSPAN,Y0）withTSPAN=[T0TFINAL]integratesthe

%systemofdifferentialequationsy'=F（t,y）fromtimeT0toTFINALwith

%initialconditionsY0.'F'isastringcontainingthenameofanODE

%file.FunctionF（T,Y）mustreturnacolumnvector.Eachrowin

%solutionarrayYcorrespondstoatimereturnedincolumnvectorT.To

%obtainsolutionsatspecifictimesT0,T1,...,TFINAL（allincreasing

%oralldecreasing）,useTSPAN=[T0T1...TFINAL].

%

%[T,Y]=ODE45（'F',TSPAN,Y0,OPTIONS）solvesasabovewithdefault

%integrationparametersreplacedbyvaluesinOPTIONS,anargument

%createdwiththeODESETfunction.SeeODESETfordetails.Commonly

%usedoptionsarescalarrelativeerrortolerance'RelTol'（1e-3by

%default）andvectorofabsoluteerrortolerances'AbsTol'（all

%components1e-6bydefault）.

%

%[T,Y]=ODE45（'F',TSPAN,Y0,OPTIONS,P1,P2,...）passestheadditional

%parametersP1,P2,...totheODEfileasF（T,Y,FLAG,P1,P2,...）（see

%ODEFILE）.UseOPTIONS=[]asaplaceholderifnooptionsareset.

%

%ItispossibletospecifyTSPAN,Y0andOPTIONSintheODEfile（see

%ODEFILE）.IfTSPANorY0isempty,thenODE45callstheODEfile

%[TSPAN,Y0,OPTIONS]=F（[],[],'init'）toobtainanyvaluesnotsupplied

%intheODE45argumentlist.Emptyargumentsattheendofthecalllist

%maybeomitted,e.g.ODE45（'F'）.

%

%ODE45cansolveproblemsM（t,y）*y'=F（t,y）withamassmatrixMthatis

%nonsingular.UseODESETtosetMassto'M','M（t）',or'M（t,y）'ifthe

%ODEfileiscodedsothatF（T,Y,'mass'）returnsaconstant,

%time-dependent,ortime-andstate-dependentmassmatrix,respectively.

%ThedefaultvalueofMassis'none'.ODE15SandODE23Tcansolveproblems

%withsingularmassmatrices.

%

%[T,Y,TE,YE,IE]=ODE45（'F',TSPAN,Y0,OPTIONS）withtheEventspropertyin

%OPTIONSsetto'on',solvesasabovewhilealsolocatingzerocrossings

%ofaneventfunctiondefinedintheODEfile.TheODEfilemustbe

%codedsothatF（T,Y,'events'）returnsappropriateinformation.See

%ODEFILEfordetails.OutputTEisacolumnvectoroftimesatwhich

%eventsoccur,rowsofYEarethecorrespondingsolutions,andindicesin

%vectorIEspecifywhicheventoccurred.

%

%SeealsoODEFILE

%ODE45isanimplementationoftheexplicitRunge-Kutta（4,5）pairof

%DormandandPrincecalledvariouslyRK5（4）7FM,DOPRI5,DP（4,5）andDP54.

%Itusesa"free"interpolantoforder4communicatedprivatelyby

%DormandandPrince.Localextrapolationisdone.

%DetailsaretobefoundinTheMATLABODESuite,L.F.Shampineand

%M.W.Reichelt,SIAMJournalonScientificComputing,18-1,1997.

%MarkW.ReicheltandLawrenceF.Shampine,6-14-94

%Copyright1984-2009TheMathWorks,Inc.

%$Revision:

5.74.4.10$$Date:

2009/04/2103:

24:

15$

solver_name='ode45';

%Checkinputs

ifnargin<4

options=[];

ifnargin<3

y0=[];

ifnargin<2

tspan=[];

ifnargin<1

error（'MATLAB:

ode45:

NotEnoughInputs',...

'Notenoughinputarguments.SeeODE45.'）;

end

end

end

end

%Stats

nsteps=0;

nfailed=0;

nfevals=0;

%Output

FcnHandlesUsed=isa（ode,'function_handle'）;

output_sol=（FcnHandlesUsed&&（nargout==1））;%sol=odeXX（...）

output_ty=（~output_sol&&（nargout>0））;%[t,y,...]=odeXX（...）

%Theremightbenooutputrequested...

sol=[];f3d=[];

ifoutput_sol

sol.solver=solver_name;

sol.extdata.odefun=ode;

sol.extdata.options=options;

sol.extdata.varargin=varargin;

end

%Handlesolverarguments

[neq,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）;

nfevals=nfevals+1;

%Handletheoutput

ifnargout>0

outputFcn=odeget（options,'OutputFcn',[],'fast'）;

else

outputFcn=odeget（options,'OutputFcn',@odeplot,'fast'）;

end

outputArgs={};

ifisempty（outputFcn）

haveOutputFcn=false;

else

haveOutputFcn=true;

outputs=odeget（options,'OutputSel',1:

neq,'fast'）;

ifisa（outputFcn,'function_handle'）

%WithMATLAB6syntaxpassadditionalinputargumentstooutputFcn.

outputArgs=varargin;

end

end

refine=max（1,odeget（options,'Refine',4,'fast'））;

ifntspan>2

outputAt='RequestedPoints';%outputonlyattspanpoints

elseifrefine<=1

outputAt='SolverSteps';%computedpoints,norefinement

else

outputAt='RefinedSteps';%computedpoints,withrefinement

S=（1:

refine-1）/refine;

end

printstats=strcmp（odeget（options,'Stats','off','fast'）,'on'）;

%Handletheeventfunction

[haveEventFcn,eventFcn,eventArgs,valt,teout,yeout,ieout]=...

odeevents（FcnHandlesUsed,odeFcn,t0,y0,options,varargin）;

%Handlethemassmatrix

[Mtype,M,Mfun]=odemass（FcnHandlesUsed,odeFcn,t0,y0,options,varargin）;

ifMtype>0%non-trivialmassmatrix

Msingular=odeget（options,'MassSingular','no','fast'）;

ifstrcmp（Msingular,'maybe'）

warning（'MATLAB:

ode45:

MassSingularAssumedNo',['ODE45assumes'...

'MassSingularis''no''.SeeODE15SorODE23T.']）;

elseifstrcmp（Msingular,'yes'）

error（'MATLAB:

ode45:

MassSingularYes',...

['MassSingularcannotbe''yes''forthissolver.SeeODE15S'...

'orODE23T.']）;

end

%IncorporatethemassmatrixintoodeFcnandodeArgs.

[odeFcn,odeArgs]=odemassexplicit（FcnHandlesUsed,Mtype,odeFcn,odeArgs,Mfun,M）;

f0