TREFETHEN SIAM J S CI OMPUT 2005 Society for Industrial and Applied Mathematics Vol 26 No 4 pp 12141233 Abstract A modi64257cation of the exponential timedi64256erencing fourthorder RungeKutta meth od for solving sti64256 nonlinear PDEs is present ID: 4569
Download Pdf The PPT/PDF document "FOURTHORDER TIMESTEPPING FOR STIFF PDEs ..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
A.-K.KASSAMANDL.N.TREFETHENThetermisknownastheintegratingfactor.InmanyapplicationswecanworkinFourierspaceandrenderdiagonal,sothatscalarsratherthanmatricesareinvolved.Dierentiating(1.5)givesNow,multiplying(1.2)bytheintegratingfactorgives thatis,ThishastheeectofamelioratingthestilinearpartofthePDE,andwecanuseatime-steppingmethodofourchoice(forexample,afourth-orderRunge Kuttafor-mula)toadvancethetransformedequation.Inpractice,onedoesntusetheequationasitiswrittenin(1.8),butratherreplacesactualtime,,withthetimestep,andincrementallyupdatestheformulafromonetimesteptothenext.Thisgreatlyimprovesthestability.Inboththesplitstepmethodandtheintegratingfactormethod,weuseafourth-orderRunge Kuttamethodforthetime-stepping.Thefourth-orderRunge Kuttaalgorithmthatweusedtoperformthetimeintegrationforthismethodwasc,t )(Fourth-orderRK)isthetimestepandisthenonlinearfunctionalontheright-handsideof(1.8).Forthesplitstepmethod,wesimplyreplacein(1.9)withfrom(1.4).Inarecentpaper[20],FornbergandDriscolldescribeacleverextensionoftheimplicit-explicitconceptdescribedabove.Inadditiontosplittingtheproblemintoalinearandanonlinearpart,theyalsosplitthelinearpart(aftertransformationtoFourierspace)intothreeregions:low,medium,andhighwavenum-bers.Theslidermethodinvolvesusingadierentnumericalschemeineachregion.Theadvantageofthismethodisthatonecancombinehigh-ordermethodsforthelowwavenumberswithhigh-stabilitymethodsforthehigherwavenumbers.Wecansummarizeoneversionofthismethodwiththefollowingtable. Low Medium|k| High|k| AB4/AB4 AB4/AM6 AB4/AM2* isthewavenumber,AB4denotesthefourth-orderAdams Bashforthformula,AM6denotesthesixth-orderAdams Moultonformula,andAM2denotesamodiedsecond-orderAdams Moultonformulaspeciedby 23 2Lun1 isthetimestep. A.-K.KASSAMANDL.N.TREFETHENd.Thisequationisexact,andthevariousorderETDschemescomefromhowoneap-proximatestheintegral.IntheirpaperCoxandMatthewsrstpresentasequenceofrecurrenceformulaethatprovidehigherandhigher-orderapproximationsofamulti-steptype.Theyproposeageneratingformulaistheorderofthescheme.Thecoecientsaregivenbytherecurrence 2gm1+1 3gm2+g0 CoxandMatthewsalsoderiveasetofETDmethodsbasedonRunge Kuttatime-stepping,whichtheycallETDRKschemes.Inthisreportweconsideronlythefourth-orderschemeofthistype,knownasETDRK4.AccordingtoCoxandMatthews,thederivationofthisschemeisnotatallobviousandrequiresasymbolicmanipulationsystem.TheCoxandMatthewsETDRK4formulaeare: an=eLun+L1(eLI)N(un), bn=eLun+L1(eLI)N(an+h/2), cn=eLan+L1(eLI)(2N(bn+h/2)N(un)), uneLhun+h2L3{[4Lh+eLh(43Lhh)2)]N(un) +2[2+2)+ 43Lh(Lh)2+eLh(4Lh)]N(cn+h)}. Unfortunately,inthisform,ETDRK4(andindeedanyoftheETDschemesoforderhigherthantwo)suersfromnumericalinstability.Tounderstandwhythisisthecase,considertheexpression Theaccuratecomputationofthisfunctionisawell-knownprobleminnumericalanalysisandisdiscussed,forexample,inthemonographbyHigham[25],aswellasthepaperbyFriesneretal.[21].Thereasonitisnotstraightforwardisthatforsmall(2.4)suersfromcancellationerror.WeillustratethisinTable2.1bycomparingthetruevalueof)tovaluescomputeddirectlyfrom(2.4)andfromvetermsofaTaylorexpansion.Forsmall,thedirectformulaisnogoodbecauseofcancellation,buttheTaylorpolynomialisexcellent.Forlarge,thedirectformulaisne,butthe A.-K.KASSAMANDL.N.TREFETHENThesolutionwehavefoundistoevaluate)viaanintegraloveracontourinthecomplexplanethatenclosesandiswellseparatedfrom0: 2f(t) becomesamatrixinsteadofascalar,thesameapproachworks,withtheterm1)becomingtheresolventmatrix( HerecanbeanycontourthatenclosestheeigenvaluesofContourintegralsofanalyticfunctions(scalarormatrix)inthecomplexplaneareeasytoevaluatebymeansofthetrapezoidrule,whichconvergesexponentially[15,16,24,60].Inpracticewetaketobeacircleandusuallyndthat32or64equallyspacedpointsaresucient.Whenisreal,wecanexploitthesymmetryandevaluateonlyinequallyspacedpointsontheupperhalfofacirclecenteredontherealaxis,thentaketherealpartoftheresult.ThescalarsoreigenvaluesofthatariseinadiscretizedPDEtypicallylieinornearthelefthalfofthecomplexplaneandmaycoverawiderange,whichgrowswiththespatialdiscretizationparameter.Fordiusiveproblemstheyareclosetothenegativerealaxis(e.g.,theKSequation),andfordispersiveproblemstheyareclosetotheimaginaryaxis(KdV).Suitablecontoursmayaccordinglyvaryfromproblemtoproblem.Ourexperienceshowsthatmanydierentchoicesworkwell,solongasoneiscarefultoensurethattheeigenvaluesareindeedenclosedby.Forsomediusiveproblems,itmightbeadvantageoustouseaparaboliccontourextendingto,takingadvantageofexponentialdecaydeepinthelefthalf-plane,butwehavenotusedthisapproachfortheproblemstreatedinthispaper.Fordiagonalproblems,wehavetheadditionalexibilityofbeingabletochooseacontourthatdependson,suchasacirclecenteredat.Inthisspecialcase,thecontourintegralreducessimplytoameanof)over,whichweapproximatetofullaccuracybyameanoverequallyspacedpointsalong(oragainjusttheupperhalfof,followedbytakingtherealpart).Forthedetailsofexactlyhowthiscontourintegralapproachcanbeimplemented,seetheMatlabcodeslistedinsections4and5.Thereisconsiderableexibilityaboutthisprocedure,andwedonotclaimthatourparticularimplementationsareoptimal,merelythattheyworkandareeasytoprogram.Fornondiagonalproblems,quiteabitofcomputationisinvolvedsay,32matrixinversesbutasthisisdonejustoncebeforethetime-steppingbegins(assumingthatthetimestepsareofaxedsize),theimpactonthetotalcomputingtimeissmall.Itwouldbegreaterforsomeproblemsinmultiplespacedimensions,butinsomecasesonecouldamelioratetheproblemwithapreliminarySchurfactorizationtobringthematrixtotriangularform.ContourintegralsandTaylorseriesarenottheonlysolutionsthathavebeenproposedforthisproblem.BothBeylkin[5]andFriesneretal.[21],forexample,useamethodthatisbasedonscalingandsquaring.Thatmethodisalsoeective,butthecontourintegralmethodappealstousbecauseofitsgreatergeneralityfordealingwitharbitraryfunctions.Todemonstratetheeectivenessofourstabilizationmethod,Table2.2considersthecomputationofthecoecient)of(2.5)(here=1)byuseoftheformula(2.5)andbyacontourintegralmethod.Herewefollowthesimplestapproach A.-K.KASSAMANDL.N.TREFETHENBurgersequationwithperiodicboundaryconditions, ,x,t=0)=exp(10sin03andthesimulationrunningto=1.KSequationwithperiodicboundaryconditions,conditions,,32],(3.3)u(x,t=0)=cos 1+sin withthesimulationrunningto=30.Allen CahnequationwithconstantDirichletboundaryconditions,conditions,1,1],(3.4)u(x,t=0)=47sin(001andthesimulationrunningto=3.Toimposetheboundaryconditionswedeneandworkwithhomogeneousboundaryconditionsinthevariable;thespatialdiscretizationisbyan80-pointChebyshevspectralmethod(seesection5).Weemphasizethattherstthreeproblems,becauseoftheperiodicboundaryconditions,canbereducedtodiagonalformbyFouriertransformation,whereasthefourthcannotbereducedthiswayandthusforcesustoworkwithmatricesratherthanscalars.OurresultsaresummarizedinFigures1 4.Therstplotineachgurecomparesaccuracyagainststepsizeandthusshouldbereasonablyindependentofmachineandimplementation.Ultimately,ofcourse,itiscomputertimethatmatters,andthisiswhatisdisplayedinthesecondplotineachgure,basedonourMatlabimplementationsonan800MHzPentium3machine.Otherimplementationsandothermachineswouldgivesomewhatdierentresults.BeforeconsideringthedierencesamongmethodsrevealedinFigures1 4,letusrsthighlightthemostgeneralpoint:ourcomputationsshowthatitisentirelypracticaltosolvethesedicultnonlinearPDEstohighaccuracybyfourth-ordertime-stepping.Mostsimulationsinthepasthavebeenoflowerorderintime,typicallysecondorder,butwebelievethatformostpurposes,afourth-ordermethodissuperior.TurningnowtothedierencesbetweenmethodsrevealedinFigures1 4,therstthingwenoteisthatthedierencesareveryconsiderable.Themethodsdierineciencybyfactorsasgreatas10orhigher.Wewerenotabletomakeeverymethodworkineverycase.Ifamethoddoesnotappearonagraphitmeansthatitdidnotsucceed,seeminglyforreasonsofnonlinearinstability(whichperhapsmighthavebeentackledbydealiasinginourspectraldiscretizations).Akeyfeaturetolookforintherstplotineachgureistherelativepositioningofthedierentmethods.Schemesthatarefurthertotherightforagivenaccuracytakefewerstepstoachievethataccuracy.Itispossiblethateachstepismorecostly,however,sojustbecauseaschemeachievesagoodaccuracyinfewstepsdoesnotmeanthatitisthemostecient.Thesecondplotintheguresgivesinsightintothesecomputationtimes.Wecanmakeafewcommentsoneachofthemethodsinvestigated. A.-K.KASSAMANDL.N.TREFETHEN Relative timestepRelative error at t = 30 Integrating factorETDRK4RK Slider 10 0 101 102 103 10 Computer time (s)Relative error at t = 30 Fig.3ResultsfortheKSequation.OurMatlabcodeislistedinsection 0 10 0 Relative timestepRelative error at t = 3 Split step 0 101 102 103 10 0 Computer time (s)Relative error at t = 3 Fig.4ResultsfortheAllen Cahnequation.Thisproblemisnondiagonalandmorechallengingthantheothers.OurMatlabcodeislistedinsectionproblemwithschemesofhigherthansecondorder.Forsecond-ordercalculations,splitstepschemesarecertainlycompetitive.methoddoeswellinallofthediagonalcases.Itisfast,accurate,andverystable.ThisisremarkablewhenwecompareitsperformancewiththatoftheIMEXschemesfromwhichitwasderived.Theonlyproblemwiththisschemeisthedicultyingeneralizingittonondiagonalcases.Finally,theintegratingfactorschemeperformswellfortheBurgersequation.ItdoesntdowellfortheKSequation,though,comingoworstofall,andwecouldntgetittoworkfortheKdVequationortheAllen Cahnequationwiththespatialdiscretizationthatweused.Thisisalsoalittlesurprising,consideringhowwidelyusedthisschemeis. A.-K.KASSAMANDL.N.TREFETHEN Relative timestepRelative error at t = 3 radius = 1e3 Fig.5LossofstabilityintheETDschemeappliedtotheBurgersequationastheradiusofthecontourshrinkstozero.ThecontourofzeroradiuscorrespondstoanETDcalculationdirectlyfromthedeningformula.Weusetheinitialcondition0)=cos(16)(1+sin(Astheequationisperiodic,wediscretizethespatialpartusingaFourierspectralmethod.TransformingtoFourierspacegives or,inthestandardformof(1.2),u,t denotesthediscreteFouriertransform.WesolvetheproblementirelyinFourierspaceanduseETDRK4time-steppingtosolveto=150.Figure6showstheresult,whichtooklessthan1secondofcomputertimeonourworkstation.Despitetheextraordinarysensitivityofthesolutionatlatertimestoperturbationsintheinitialdata(suchperturbationsareampliedbyasmuchas10upto=150),wearecondentthatthisimageiscorrecttoplottingaccuracy.Itwouldnothavebeenpracticaltoachievethiswithatime-steppingschemeoflowerWeproducedFigure6withtheMatlabcodelistedinFigure7.5.Anondiagonalexample:Allen Cahn.TheAllen Cahnequationisan-otherwell-knownequationfromtheareaofreaction-diusionsystems: A.-K.KASSAMANDL.N.TREFETHEN%kursiv.m-solutionofKuramoto-SivashinskyequationbyETDRK4scheme%u_t=-u*u_x-u_xx-u_xxxx,periodicBCson[0,32*pi]%computationisbasedonv=fft(u),solineartermisdiagonal%comparep27.minTrefethen,"SpectralMethodsinMATLAB",SIAM2000%AKKassamandLNTrefethen,July2002%Spatialgridandinitialcondition:N=128;x=32*pi*(1:N)/N;u=cos(x/16).*(1+sin(x/16));v=fft(u);%PrecomputevariousETDRK4scalarquantities:h=1/4;%timestepk=[0:N/2-10-N/2+1:-1]/16;%wavenumbersL=k.^2-k.^4;%FouriermultipliersE=exp(h*L);E2=exp(h*L/2);M=16;%no.ofpointsforcomplexmeansr=exp(1i*pi*((1:M)-.5)/M);%rootsofunityLR=h*L(:,ones(M,1))+r(ones(N,1),:);Q=h*real(mean((exp(LR/2)-1)./LR,2));f1=h*real(mean((-4-LR+exp(LR).*(4-3*LR+LR.^2))./LR.^3,2));f2=h*real(mean((2+LR+exp(LR).*(-2+LR))./LR.^3,2));f3=h*real(mean((-4-3*LR-LR.^2+exp(LR).*(4-LR))./LR.^3,2));%Maintime-steppingloop:uu=u;tt=0;tmax=150;nmax=round(tmax/h);nplt=floor((tmax/100)/h);g=-0.5i*k;forn=1:nmaxt=n*h;Nv=g.*fft(real(ifft(v)).^2);a=E2.*v+Q.*Nv;Na=g.*fft(real(ifft(a)).^2);b=E2.*v+Q.*Na;Nb=g.*fft(real(ifft(b)).^2);c=E2.*a+Q.*(2*Nb-Nv);Nc=g.*fft(real(ifft(c)).^2);v=E.*v+Nv.*f1+2*(Na+Nb).*f2+Nc.*f3;ifmod(n,nplt)==0u=real(ifft(v));uu=[uu,u];tt=[tt,t];%Plotresults:surf(tt,x,uu),shadinginterp,lightingphong,axistightview([-9090]),colormap(autumn);set(gca,zlim,[-550])light(color,[110],position,[-1,2,2])material([0.300.600.6040.001.00]);Fig.7MatlabcodetosolvetheKSequationandproduceFigure.Despitetheextraordinarysensitivityofthisequationtoperturbations,thiscodecomputescorrectresultsinlessthansecondonan800MHzPentiummachine. A.-K.KASSAMANDL.N.TREFETHEN%allencahn.m-solutionofAllen-CahnequationbyETDRK4scheme%u_t=0.01*u_xx+u-u^3on[-1,1],u(-1)=-1,u(1)=1%computationisbasedonChebyshevpoints,solineartermisnondiagonal%comparep34.minTrefethen,"SpectralMethodsinMATLAB",SIAM2000%AKKassamandLNTrefethen,July2002%Spatialgridandinitialcondition:N=20;[D,x]=cheb(N);x=x(2:N);%spectraldifferentiationmatrixw=.53*x+.47*sin(-1.5*pi*x)-x;%usew=u-xtomakeBCshomogeneousu=[1;w+x;-1];%PrecomputevariousETDRK4matrixquantities:h=1/4;%timestepM=32;%no.ofpointsforresolventintegralr=15*exp(1i*pi*((1:M)-.5)/M);%pointsalongcomplexcircleL=D^2;L=.01*L(2:N,2:N);%2nd-orderdifferentiationA=h*L;E=expm(A);E2=expm(A/2);I=eye(N-1);Z=zeros(N-1);f1=Z;f2=Z;f3=Z;Q=Z;forj=1:Mz=r(j);zIA=inv(z*I-A);Q=Q+h*zIA*(exp(z/2)-1);f1=f1+h*zIA*(-4-z+exp(z)*(4-3*z+z^2))/z^2;f2=f2+h*zIA*(2+z+exp(z)*(z-2))/z^2;f3=f3+h*zIA*(-4-3*z-z^2+exp(z)*(4-z))/z^2;f1=real(f1/M);f2=real(f2/M);f3=real(f3/M);Q=real(Q/M);%Maintime-steppingloop:uu=u;tt=0;tmax=70;nmax=round(tmax/h);nplt=floor((tmax/70)/h);forn=1:nmaxt=n*h;Nu=(w+x)-(w+x).^3;a=E2*w+Q*Nu;Na=a+x-b=E2*w+Q*Na;Nb=b+x-c=E2*a+Q*(2*Nb-Nu);Nc=c+x-w=E*w+f1*Nu+2*f2*(Na+Nb)+f3*Nc;ifmod(n,nplt)==0u=[1;w+x;-1];uu=[uu,u];tt=[tt,t];%Plotresults:surf([1;x;-1],tt,uu),lightingphong,axistightview([-4560]),colormap(cool),light(col,[110],pos,[-10010])Fig.9MatlabcodetosolvetheAllen CahnequationandproduceFigure.Again,thiscodetakeslessthansecondtorunonan800MHzPentiummachine. A.-K.KASSAMANDL.N.TREFETHENTREFETHENA.IserlesAFirstCourseintheNumericalAnalysisofDierentialEquations,CambridgeUniversityPress,Cambridge,UK,2000.2000.A.IserlesOnthenumericalquadratureofhighlyoscillatingintegrals:FouriertransformsIMAJ.Numer.Anal.,24(2004),pp.365 391.365 391.J.C.Jiminez,R.Biscay,C.Mora,andL.M.RodriguezDynamicpropertiesofthelocallin-earizationmethodforinitial-valueproblems,Appl.Math.Comput.,126(2002),pp.63 81.63 81.G.E.Karniadakis,M.Israeli,andS.A.OrszagHighordersplittingmethodsfortheincompressibleNavier-Stokesequations,J.Comput.Phys.,97(1991),pp.414 443.414 443.C.A.KennedyandM.H.CarpenterAdditiveRunge Kuttaschemesforconvection-diusion-reactionequations,Appl.Numer.Math.,44(2003),pp.139 181.139 181.J.KimandP.MoinApplicationsofafractionalstepmethodtoincompressibleNavier-Stokesequations,J.Comput.Phys.,59(1985),pp.308 323.308 323.O.M.Knio,H.N.Najim,andP.S.WyckoffAsemi-implicitnumericalschemeforreacting,J.Comput.Phys.,154(1999),pp.428 467.428 467.D.KortewegandG.deVriesOnthechangeofformoflongwavesadvancinginarectan-gularcanal,andonanewtypeoflongstationarywaves,Philos.Mag.Ser.5,39(1895),pp.422 433.422 433.W.KressandB.GustafssonDeferredcorrectionmethodsforinitialvalueboundaryprob-,J.Sci.Comput.,17(2002),pp.241 251.241 251.M.KrusemeyerDierentialEquations,MacmillanCollegePublishing,NewYork,1994.1994.Y.KuramotoandT.TsuzukiPersistentpropagationofconcentrationwavesindissipativemediafarfromthermalequilibrium,Prog.Theoret.Phys.,55(1976),pp.356 369.356 369.Y.Maday,A.T.Patera,andE.M.RønquistAnoperator-integration-factorsplittingmethodfortime-dependentproblems:Applicationtoincompressibleuidow,J.Sci.Com-put.,5(1990),pp.263 292.263 292.R.I.McLachlanandP.AtelaTheaccuracyofsymplecticintegrators,Nonlinearity,5(1992),pp.541 562.541 562.R.McLachlanSymplecticintegrationofHamiltonianwaveequations,Numer.Math.,66(1994),pp.465 492.465 492.W.J.MerryfieldandB.ShizgalPropertiesofcollocationthird-derivativeoperators,J.Comput.Phys.,105(1993),pp.182 185.182 185.P.A.MilewskiandE.G.TabakApseudospectralprocedureforthesolutionofnonlinearwaveequationswithexamplesfromfree-surfaceows,SIAMJ.Sci.Comput.,21(1999),pp.1102 1114.1102 1114.M.L.MinionSemi-implicitspectraldeferredcorrectionmethodsforordinarydierentialequa-,Commun.Math.Sci.,1(2003),pp.471 500.471 500.D.R.Mott,E.S.Oran,andB.vanLeerAquasi-steadystatesolverforthestiordinarydierentialequationsofreactionkinetics,J.Comput.Phys.,164(2000),pp.407 428.407 428.B.Nicolaenko,B.Scheurer,andT.TemamSomeglobalpropertiesoftheKuramoto-Sivashinskyequation:Nonlinearstabilityandattractors,Phys.D,16(1985),pp.155 183.155 183.S.P.NørsettAnA-stablemodicationoftheAdams Bashforthmethods,inConferenceonNumericalSolutionofDierentialEquations(Dundee,1969),LectureNotesinMath.109,Springer-Verlag,Berlin,1969,pp.214 219.214 219.E.OttChaosinDynamicalSystems,CambridgeUniversityPress,Cambridge,UK,1993.1993.R.D.RuthAcanonicalintegrationtechnique,IEEETrans.NuclearScience,NS-30(1983),pp.2669 2671.2669 2671.S.J.RuuthImplicit-explicitmethodsforreaction-diusionproblemsinpatternformation,J.Math.Biol.,34(2)(1995),pp.148 176.148 176.Y.SaadAnalysisofsomeKrylovsubspaceapproximationstothematrixexponentialoperatorSIAMJ.Numer.Anal.,29(1992),pp.209 228.209 228.J.M.Sanz-SernaandM.P.CalvoNumericalHamiltonianProblems,ChapmanandHall,London,1994.1994.M.SchatzmanTowardnon-commutativenumericalanalysis:HighorderintegrationintimeJ.Sci.Comput.,17(2002),pp.99 116.99 116.L.M.SmithandF.WaleffeTransferofenergytotwo-dimensionallargescalesinforced,rotatingthree-dimensionalturbulence,Phys.Fluids,11(1999),pp.1608 1622.1608 1622.L.M.SmithandF.WaleffeGenerationofslowlargescalesinforcedrotatingstratiedturbulence,J.FluidMech.,451(2002),pp.145 169.145 169.G.StrangOntheconstructionandcomparisonofdierenceschemes,SIAMJ.Numer.Anal.,5(1968),pp.506 517.506 517.E.TadmorTheexponentialaccuracyofFourierandChebyshevdierencingmethods,SIAMJ.Numer.Anal.,23(1986),pp.1 10.