A brief introduction to using ode45 in MATLAB MATLAB中使用ode45的简单介绍

PAGE

1

NurAdilaFarukSenan

DepartmentofMechanicalEngineeringUniversityofCaliforniaatBerkeley

Abriefintroductiontousingode45inMATLAB

MATLAB’sstandardsolverforordinarydifferentialequations(ODEs)isthefunctionode45.ThisfunctionimplementsaRunge-Kuttamethodwithavariabletimestepforefficientcomputation.ode45isdesignedtohandlethefollowinggeneralproblem:

dx=f(t,x), x(t)=x, (1)

dt 0 0

wheretistheindependentvariable,xisavectorofdependentvariablestobefoundandf(t,x)isafunctionoftandx.Themathematicalproblemisspecifiedwhenthevectoroffunctionsontheright-handsideofEq.(1),f(t,x),issetandtheinitialconditions,x=x0attimet0aregiven.

InME175,thesolutionisoftennotcompleteonceyouhavesolvedtheproblemandobtainedtheode’sgoverningthesystemsmotion.Itisoftenbeneficialtoproduceavisualrepresentationofwhatexactlythetrajectoriesofaparticlerepresentedbyahighlycomplicatedlookingordinarydifferentialequationlookslikeandthefollowingisabriefexplanationofhowtoaccomplishthis.

Step1:Reducethegivenodetoaseriesoffirstorderequations

Thisistheveryfirststepyoushoulddo,ideallyonaseparatepieceofpaper.Forexample,ifthegivenodeis,

my¨+y˙ey−y2=5, y0=3,y˙0=−1, (2)thenthevectorxhastwocomponents:yandy˙,or,

x(1)=y

and

x(2)=

dx(1)=x(2)dt

y˙, (3)

x(2)= (5−x(2)ex(1)+(x(1))2), (4)

d 1

dt

m

wherewehavemadeuseoftherepresentationsfory,y˙andy¨giveninEqn(3)towrite(4).

Whatifwehavemorethanoneode?Forexample,supposethatinadditiontoEq.(2)above,wealsohaveasecondequation,

d3z+d2z

dt3

dt2−sin(z)=t, z0=0,z˙0=0,z¨0=1. (5)

Then,wecaneasilysolveforbothyandzsimultaneouslybymakingxalargervector:

x(3)=z

dz

x(4)=

x(5)=

dtd2zdt2

(6)

Thus,

and

x=[y,y˙,z,z˙,z¨], (7)

|t=0

y|t=0,dt|t=0,z|t=0,dt|t=0,dt2|t=0

x =＿ dy

dz d2z ＿

(8)

.

Thereisnopreferenceforplacingthey-relatedvariablesasthefirsttwovariablesinsteadofthez-relatedvariables.Infact,anyarbitraryorderisperfectlyvalid,suchas,

x=[y,z,y˙,z˙,z¨] and x=[z,z˙,z¨,y,y˙]. (9)

Whatisimportanthoweveristhatyouremainconsistentthroughoutyourcomputationsintermsofhowyouwriteoutyoursystemoffirstorderode’s.1Thus,ifyouusetheorderinggiveninEq.(7),thenthesystemoffirstorderode’sarecomprisedof,

d

x(1)=x(2)

dt

d 1( )

x(2)= 5−x(2)ex(1)+(x(1))2,

dt m

dx(3)=x(4)dt

dx(4)=x(5)dt

dx(5)=t−x(5)+sin(x(3)), (10)

dt

whileemployingtherepresentationx=[y,z,y˙,z˙,z¨]resultsin,

d

x(1)=x(3)

dt

dx(2)=x(4)dt

d 1( )

x(3)= 5−x(2)ex(2)+(x(2))2,

dt m

dx(4)=x(5)dt

dx(5)=t−x(5)+sin(x(2)). (11)

dt

Basically,youcouldhaveanarbitrarynumberofhigherorderode’s.Whatisimportantisthatyoureducethemtomultiplefirstorderode’sandkeeptrackofwhatorderyou’veassignedthedifferentderivativesinthefinalxvectorthatwillbesolved.

Step2:Codingitin

Nowthatyouhaveeverythinginfirstorderform,youwillneedthefollowingcommandsinyourmaincode:

[t,x]=ode45(@fname,tspan,xinit,options)

fnameisthenameofthefunctionMfileusedtoevaluatetheright-hand-sidefunctioninEq.(1).Thisisthefunctionwherewewillinputthesystemoffirstorderode’stobeintegrated(suchasinEqs.(10)and(11)).Iwillexplainthisinalittlemoredetaillateron.

1Ofcourse,theremightbesomesubtletieswithregardstohowode45numericallyintegratesthegivenequationandinsomecases,itmightmakesensetosolvesomeequationsbeforeothersbutforthesimpleproblemswewillbedealingwith,thisisanon-issue.

tspanisthevectordefiningthebeginningandendlimitsofintegration,aswellashowlargewewantourtimestepstobe.I.e.ifweareintegratingfromt=0tot=10,andwanttotake100timesteps,thentspan=[0:0.1:10]ortspan=linspace(0,10,100)).

xinitisthevectorofinitialconditions.Makesurethattheordercorrespondstotheorderingusedtowritey,zandtheirderivativesintermsofx.Alsonotethatifxconsistsof5variables,thenweneedaninputof5initialconditions(seeEqn.(8)).

optionsissomethingthatisverywellexplainedinthehelpsessionofMATLAB.Formostpurposes,usingthedefaultvalueissufficient.

tisthevalueoftheindependentvariableatwhichthesolutionarrayxiscal-culated.Thisvectorisnotnecessarilyequaltotspanabovebecauseode45doessomeamountofplayingaboutwithstepsizestomaximizebothaccuracyandef-ficiency(takingsmallerstepswheretheproblemchangesrapidlyandlargerstepswhereit’srelativelyconstant).Thelengthofthoweveristhesameasthatoftspan

xiswhatthisisallabout.2xisanarray(ormatrix)withsizelength(t)bylength(xinit).Eachcolumnofxisadifferentdependentvariable.Forexample,ifwehavex=[y,y˙,z,z˙,z¨],andassume(forsimplicity)thatwerequireonlythevaluesatt=0,1,2,3,...,10(i.e.weevaluatethefunctionat11differenttimes)then

x=

y|t=0

y|t=1

y|t=10

y˙|t=0 z|t=0

y˙|t=1 z|t=1

...

...

...

y˙|t=10 z|t=10

z˙|t=0

z˙|t=1

z˙|t=10

z¨|t=0z¨|t=1

z¨|t=10

. (12)

Thus,x(1,4)givesusthevalueofz˙att=0,x(7,4)givesusthevalueofz˙att=6,andx(11,4)givesusthevalueofz˙att=10sincez˙isthefourthvariableinx.Thesameholdsforalltheothervariables.Inshort,

x(:,k)givesusthecolumnvectorofthek’thvariable:k=1wouldcorrespondto

y,k=2toy˙,etc.

x(j,:)givesusthevaluesofallthecomputedvariablesatasingleinstantoftimecorrespondingtoindexj.

2BeforetheinventionoftheHokeyPokeydance,ancientchildrenattendingancientcampfireeventswouldsitaroundancientcampfiressinging:

YouputyourleftfootinYouputyourleftfootoutYouputyourleftfootinAndyoushakeitallabout

Youusetheode45MATLABfunctiontogetyourhomeworksdoneontimexiswhatit’sallabout!

Unfortunately,thelackofcomputers(andMATLAB)soonmadethisversionobsolete.

⇒x(j,k)≡x(time,variableofinterest).

Thelinefollowing[t,x]=ode45(@fname,tspan,xinit,options)shouldbeare-definitionofyourvariables.Thus,ifyouusedx=[y,y˙,z,z˙,z¨],thenyouwouldwrite:

[t,x]=ode45(fname,tspan,xinit,options)y=x(:,1);

ydot=x(:,2);

z=x(:,3);

zdot=x(:,4);

zdotdot=x(:,5);

Ofcourse,youcouldjustaswellnotbothertodefiney,ydot,etc.andinsteaddealdirectlywiththeindicialformofx(e.g.usingx(:,1)wheneverwemeany,etc)butdefiningthevariablesclearlydoesmakethingsalittleeasiertodebugat3inthemorningthedaythehomeworkisdue.

Belowthis,youusuallyplotwhatexactlyitisyou’reinterestedinshowing(orhavebeenaskedtoshow):Thetrajectoryofaparticleasafunctionoftime?Therelationshipbetweenrandθforaplanarpendulum?etc.TheplotandsubplotcommandsinMATLABarelucidlyexplainedintheMATLABhelpandIwon’tgointodetailaboutthemhere.Bearinmindthatifyouplantohandin20plots,youwilldothegrader(andmothernature)afavorbyusingthesubplotfunctiontofitmultipleplotsintoonepage.Don’tgooverboardwiththishowever-20plotsonasinglepageisn’tagoodideaeither.Additionally,don’tforgettolabelyourgraphsadequately:Eachplotshouldhaveatitle,labelsforthex-axis,y-axisandalegendiftherearemultiplecurvesonthesameplot.Aplotwithoutattachedlablesisjustabunchoflines!

Finally,notethateverythingin[t,x]=ode45(@fname,tspan,xinit,options)arejustvariables.Youcanusewhatevertermsyouwish-Tinsteadoft,x0insteadofxinit,orOUTinsteadofxtociteafewexamples.Justrememberthatifyouhavedefinedavariable,thenyouhavetoalwaysrefertoitbythesamename.

[TIME,OUT]=ode45(@...)followedbyy=x(:,1);willgiveyouanerrorsincexisundefined.(Thecorrectequationshouldbey=OUT(:,1);)

Anothercommonmistakeistoredefinethesamevariable.Forexample,doingsome-thinglikex=x(:,5);.Ifyou’relucky,youwillgetanerrormessage.Ifyou’renothowever(whichusuallyhappensiftheerrorisnotasglaring),thenyoucouldendupspendinghourstryingtofigureoutwhyyoursolutionlooksalittlefunny.

Aside:Whatis’fname?

Recallthatwestillhaven’ttoldMATLABwhatexactlytheequationsofmotionarethatneedtobeintegrated.Thisiswherefnamecomesin.fnameisthenameofthefunctioncontainingallthefirstorderode’swewroterightatthebeginning.Wecannamethisfunctionanythingwelikesolongasthenameyougiveitisthesameaswhatyouusewhencallingitin[t,x]=ode45(@fname,tspan,xinit,options).

Thus,ifyoucallyourfunctionsuperman,then

[t,x]=ode45(@superman,tspan,xinit,options)isrightbut

[t,x]=ode45(@batman,tspan,xinit,options)iswrong(ifyouhaveafunction

namedbatmaninyourlibrary,thenyou’llendupgettingthesolutiontothatprobleminstead!).

Furthermore,,thefunctiondoesn’thavetobeinthesamemfileasyouroriginalcode

-somepeopleprefertowriteitasasub-functionrightattheendoftheprogram,especiallyifthecodeisn’ttoolargeorcomplicated.Forexample,sayyourfunctioniscalledME175example.Then,yourmfilewouldlooklikethis:

functiondxdt=ME175example(t,x)

%Here,t,xanddxdtareagain,justvariables.Youcancallthemwhateveryouwant,

%solongastistheindependentvariable,xisthevectorofdependentvariablesand

%dxdtisthevectoroffirstorderderivatives(withrespecttotheindependentvariable)

%Definetheconstants:

m=1;

%Definethevariablestomakethingslittlemorelegible

%Recallthatx=[y,ydot,z,zdot,zdotdot]

y=x(1);

ydot=x(2);z=x(3);

zdot=x(4);zdotdot=x(5);

%Noticethathere,xisjusta1by5rowvectorandnotalargearraylikein

%[t,x]=ode45(@fname,tspan,xinit,options)

%thearraydxdtisthesamelengthasx

dxdt=zeros(size(x));

( ）

dxdt(1)=ydot;

dxdt(2)=1/m5−x(2)exp(y)+y2; %Thisisydotdot

dxdt(3)=zdot;

dxdt(4)=zdotdot;

dxdt(5)=t-zdotdot+sin(z); %Thisiszdotdotdot

%Notethattheinputargumentsmustbetandx(inthatorder)eveninthecasewheretisnotexplicitlyusedinthefunction.

Abasictemplate

ThefollowingisabasictemplateyoucanjustcutandpasteintoMATLABeverytimeyouwouldliketointegrateahigherorderODE:

functioncallthiswhateveryouwant

%definetstart,tendandthenumberoftimestepsyouwanttotake,n:tstart=?;

tend=?;n=?;

tspan=linspace(tstart,tend,n);

%definetheinitialconditionsmakingsuretousetherightorderingxinit=[...;...;...;...];

%Getx.Ihaveleftoptionsasthedefaultvaluesoit’snotincludedinthefunctioncallbelow[t,x]=ode45(@integratingfunction,tspan,xinit)

%definetheoutputvariables:

y=x(:,1);

ydot=x(:,2);etc...

%Includewhateverplottingfunctionsyouneedtoincludehere.Youcanplotanythingyouwant:%y(t),y˙(t),y˙(y),y˙(y¨),etc.

subplot(?,?,?)

plot(whateveryouareplotting)

%plot3(...)%usethisisyouare