6 ORDINARY DIFFERENTIAL EQUATIONS.docx

6 ORDINARY DIFFERENTIAL EQUATIONS.docx

《6 ORDINARY DIFFERENTIAL EQUATIONS.docx》由会员分享，可在线阅读，更多相关《6 ORDINARY DIFFERENTIAL EQUATIONS.docx（112页珍藏版）》请在冰豆网上搜索。

6ORDINARYDIFFERENTIALEQUATIONS

6ORDINARYDIFFERENTIALEQUATIONS

Differentialequationsaremathematicaldescriptionsofhowthevariablesandtheirderivatives（ratesofchange）withrespecttooneormoreindependentvariableaffecteachotherinadynamicalway.Theirsolutionsshowushowthedependentvariable（s）willchangewiththeindependentvariable（s）.Manyproblemsinnaturalsciencesandengineeringfieldsareformulatedintoascalardifferentialequationoravectordifferentialequation—thatis,asystemofdifferentialequations.

Inthischapter,welookintoseveralmethodsofobtainingthenumericalsolutionstoordinarydifferentialequations（ODEs）inwhichalldependentvariables（x）dependonasingleindependentvariable（t）.First,theinitialvalueproblems（IVPs）willbehandledwithseveralmethodsincludingRunge–Kuttamethodandpredictor–correctormethodsinSections6.1to6.5.Thefinalsection（Section6.6）willintroducetheshootingmethodandthefinitedifferenceandcollocationmethodforsolvingthetwo-pointboundaryvalueproblem（BVP）.ODEsarecalledanIVPifthevaluesx（t0）ofdependentvariablesaregivenattheinitialpointt0oftheindependentvariable,whiletheyarecalledaBVPifthevaluesx（t0）andx（tf）aregivenattheinitial/finalpointst0andtf.

6.1EULER’SMETHOD

WhentalkingaboutthenumericalsolutionstoODEs,everyonestartswiththeEuler’smethod,sinceitiseasytounderstandandsimpletoprogram.EventhoughitslowaccuracykeepsitfrombeingwidelyusedforsolvingODEs,itgivesusacluetothebasicconceptofnumericalsolutionforadifferentialequationsimplyandclearly.Let’sconsiderafirst-orderdifferentialequation:

Ithasthefollowingformofanalyticalsolution:

whichcanbeobtainedbyusingaconventionalmethodortheLaplacetransformtechnique.However,suchaniceanalyticalsolutiondoesnotexistforeverydifferentialequation;evenifitexists,itisnoteasytofindevenbyusingacomputerequippedwiththecapabilityofsymboliccomputation.Thatiswhyweshouldstudythenumericalsolutionstodifferentialequations.

Then,howdowetranslatethedifferentialequationintoaformthatcaneasilybehandledbycomputer?

Firstofall,wehavetoreplacethederivativey’（t）=dy/dtinthedifferentialequationbyanumericalderivative（introducedbefore）,wherethestep-sizehisdeterminedbasedontheaccuracyrequirementsandthecomputationtimeconstraints.Euler’smethodapproximatesthederivativeinEq.（6.1.1）withEq.（5.1.2）as

andsolvesthisdifferenceequationstep-by-stepwithincreasingtbyheachtimefromt=0.

Thisisanumericsequence{y（kh）},whichwecallanumericalsolutionofEq.（6.1.1）.

Tobespecific,lettheparametersandtheinitialvalueofEq.（6.1.1）bea=1,r=1,andy0=0.Then,theanalyticalsolution（6.1.2）becomes

y（t）=1−e−at（6.1.5）

andthenumericalsolution（6.1.4）withthestep-sizeh=0.5andh=0.25areaslistedinTable6.1anddepictedinFig.6.1.WemakeaMATLABprogram“nm610.m”,whichusesEuler’smethodforthedifferentialequation（6.1.1）,actuallysolvingthedifferenceequation（6.1.3）andplotsthegraphsofthenumericalsolutionsinFig.6.1.

%nm610:

Eulermethodtosolvea1st-orderdifferentialequation

clear,clf

a=1;r=1;y0=0;tf=2;

t=[0:

0.01:

tf];yt=1-exp（-a*t）;%Eq.（6.1.5）:

trueanalyticalsolution

plot（t,yt,’k’）,holdon

klasts=[842];hs=tf./klasts;

y

（1）=y0;

foritr=1:

3%withvariousstepsizeh=1/8,1/4,1/2

klast=klasts（itr）;h=hs（itr）;y

（1）=y0;

fork=1:

klast

y（k+1）=（1-a*h）*y（k）+h*r;%Eq.（6.1.3）:

plot（[k-1k]*h,[y（k）y（k+1）],’b’,k*h,y（k+1）,’ro’）

ifk<4,pause;end

end

end

Figure6.1ExamplesofnumericalsolutionobtainedbyusingtheEuler’smethod.

Thegraphsseemtotellusthatasmallstep-sizehelpsreducetheerrorsoastomakethenumericalsolutionclosertothe（true）analyticalsolution.But,aswillbeinvestigatedthoroughlyinSection6.2,itisonlypartiallytrue.Infact,atoosmallstep-sizenotonlymakesthecomputationtimelonger（proportionalas1/h）,butalsoresultsinratherlargererrorsduetotheaccumulatedround-offeffect.Thisiswhyweshouldlookforothermethodstodecreasetheerrorsratherthansimplyreducethestep-size.

Euler’smethodcanalsobeappliedforsolvingafirst-ordervectordifferentialequation

whichisequivalenttoahigh-orderscalardifferentialequation.Thealgorithmcanbedescribedby

andiscastintotheMATLABroutine“ode_Euler（）”.

function[t,y]=ode_Euler（f,tspan,y0,N）

%Euler’smethodtosolvevectordifferentialequationy’（t）=f（t,y（t））

%fortspan=[t0,tf]andwiththeinitialvaluey0andNtimesteps

ifnargin<4|N<=0,N=100;end

ifnargin<3,y0=0;end

h=（tspan

（2）-tspan

（1））/N;%stepsize

t=tspan

（1）+[0:

N]’*h;%timevector

y（1,:

）=y0（:

）’;%alwaysmaketheinitialvaluearowvector

fork=1:

N

y（k+1,:

）=y（k,:

）+h*feval（f,t（k）,y（k,:

））;%Eq.（6.1.7）

end

6.2HEUN’SMETHOD:

TRAPEZOIDALMETHOD

Anothermethodofsolvingafirst-ordervectordifferentialequationlikeEq.（6.1.6）comesfromintegratingbothsidesoftheequation.

Ifweassumethatthevalueofthe（derivative）functionf（t,y）isconstantasf（tk,y（tk））withinonetimestep[tk,tk+1）,thisbecomesEq.（6.1.7）（withh=tk+1−tk）,amountingtoEuler’smethod.Ifweusethetrapezoidalrule（5.5.3）,itbecomes

function[t,y]=ode_Heun（f,tspan,y0,N）

%Heunmethodtosolvevectordifferentialequationy’（t）=f（t,y（t））

%fortspan=[t0,tf]andwiththeinitialvaluey0andNtimesteps

ifnargin<4|N<=0,N=100;end

ifnargin<3,y0=0;end

h=（tspan

（2）-tspan

（1））/N;%stepsize

t=tspan

（1）+[0:

N]’*h;%timevector

y（1,:

）=y0（:

）’;%alwaysmaketheinitialvaluearowvector

fork=1:

N

fk=feval（f,t（k）,y（k,:

））;y（k+1,:

）=y（k,:

）+h*fk;%Eq.（6.2.3）

y（k+1,:

）=y（k,:

）+h/2*（fk+feval（f,t（k+1）,y（k+1,:

）））;%Eq.（6.2.4）

end

But,theright-handside（RHS）ofthisequationhasyk+1,whichisunknownattk.Toresolvethisproblem,wereplacetheyk+1ontheRHSbythefollowingEuler’sapproximation:

sothatitbecomes

ThisisHeun’smethod,whichisimplementedintheMATLABroutine“ode_Heun（）”.Itisakindofpredictor-and-correctormethodinthatitpredictsthevalueofyk+1byEq.（6.2.3）attkandthencorrectsthepredictedvaluebyEq.（6.2.4）attk+1.ThetruncationerrorofHeun’smethodisO（h2）（proportionaltoh2）asshowninEq.（5.6.1）,whiletheerrorofEuler’smethodisO（h）.

6.3RUNGE–KUTTAMETHOD

AlthoughHeun’smethodisalittlebetterthantheEuler’smethod,itisstillnotaccurateenoughformostreal-worldproblems.Thefourth-orderRunge–Kutta（RK4）methodhavingatruncationerrorofO（h4）isoneofthemostwidelyusedmethodsforsolvingdifferentialequations.Itsalgorithmisdescribedbelow.

where

Equation（6.3.1）isthecoreofRK4method,whichmaybeobtainedbysubstitutingSimpson’srule（5.5.4）

intotheintegralform（6.2.1）ofdifferentialequationandreplacingfk+1/2withtheaverageofthesuccessivefunctionvalues（fk2+fk3）/2.Accordingly,theRK4methodhasatruncationerrorofO（h4）asEq.（5.6.2）andthusisexpectedtoworkbetterthantheprevioustwomethods.

Thefourth-orderRunge–Kutta（RK4）methodiscastintotheMATLABroutine“ode_RK4（）”.Theprogram“nm630.m”usesthisroutinetosolveEq.（6.1.1）withthestepsizeh=（tf−t0）/N=2/4=0.5andplotsthenumericalresulttogetherwiththe（true）analyticalsolution.

function[t,y]=ode_RK4（f,tspan,y0,N,varargin）

%Runge-Kuttamethodtosolvevectordifferentialeqny’（t）=f（t,y（t））

%fortspan=[t0,tf]andwiththeinitialvaluey0andNtimesteps

ifnargin<4|N<=0,N=100;end

ifnargin<3,y0=0;end

y（1,:

）=y0（:

）’;%makeitarowvector

h=（tspan

（2）-tspan

（1））/N;t=tspan

（1）+[0:

N]’*h;

fork=1:

N

f1=h*feval（f,t（k）,y（k,:

）,varargin{:

}）;f1=f1（:

）’;%（6.3.2a）

f2=h*feval（f,t（k）+h/2,y（k,:

）+f1/2,varargin{:

}）;f2=f2（:

）’;%（6.3.2b）

f3=h*feval（f,t（k）+h/2,y（k,:

）+f2/2,varargin{:

}）;f3=f3（:

）’;%（6.3.2c）

f4=h*feval（f,t（k）+h,y（k,:

）+f3,varargin{:

}）;f4=f4（:

）’;%（6.3.2d）

y（k+1,:

）=y（k,:

）+（f1+2*（f2+f3）+f4）/6;%Eq.（6.3.1）

end

%nm630:

Heun/Euer/RK4methodtosolveadifferentialequation（d.e.）

clear,clf

tspan=[02];

t=tspan

（1）+[0:

100]*（tspan

（2）-tspan

（1））/100;

a=1;yt=1-exp（-a*t）;%Eq.（6.1.5）:

trueanalyticalsolution

plot（t,yt,’k’）,holdon

df61=inline（’-y+1’,’t’,’y’）;%Eq.（6.1.1）:

d.e.tobesolved

y0=0;N=4;

[t1,ye]=oed_Euler（df61,tspan,y0,N）;

[t1,yh]=ode_Heun（df61,tspan,y0,N）;

[t1,yr]=ode_RK4（df61,tspan,y0,N）;

plot（t,yt,’k’,t1,ye,’b:

’,t1,yh,’b:

’,t1,yr,’r:

’）

plot（t1,ye,’bo’,t1,yh,’b+’,t1,yr,’r*’）

N=1e3;%toestimatethetimeforNiterations

tic,[t1,ye]=ode_Euler（df61,tspan,y0,N）;time_Euler=toc

tic,[t1,yh]=ode_Heun（df61,tspan,y0,N）;time_Heun=toc

tic,[t1,yr]=ode_RK4（df61,tspan,y0,N）;time_RK4=toc

ComparisonofthisresultwiththoseofEuler’smethod（“ode_Euler（）”）andHeun’smethod（“ode_Heun（）”）isgiveninFig.6.2,whichshowsthattheRK4methodisbetterthanHeun’smethod,whileEuler’smethodistheworstintermsofaccuracywiththesamestep-size.But,

Figure6.2Numericalsolutionsforafirst-orderdifferentialequation.

intermsofcomputationalload,theorderisreversed,becauseEuler’smethod,Heun’smethod,andtheRK4methodneed1,2,and4functionevaluations（calls）periteration,respectively.

（cf）Notethatafunctioncalltakesmuchmoretimethanamultiplicationandthusthenumberoffunctioncallsshouldbeacriterioninestimatingandcomparingcomputationaltime.

TheMATLABbuilt-inroutines“ode23（）”and“ode45（）”implementtheRunge–Kutt