106K - views

FOURTHORDER TIMESTEPPING FOR STIFF PDEs ALYKHAN KASSAM AND LLOYD N

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

Tags : TREFETHEN SIAM
Download Pdf

FOURTHORDER TIMESTEPPING FOR STIFF PDEs ALYKHAN KASSAM AND LLOYD N




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.



Presentation on theme: "FOURTHORDER TIMESTEPPING FOR STIFF PDEs ALYKHAN KASSAM AND LLOYD N"— Presentation transcript:

A.-K.KASSAMANDL.N.TREFETHENThetermisknownastheintegratingfactor.InmanyapplicationswecanworkinFourierspaceandrenderdiagonal,sothatscalarsratherthanmatricesareinvolved.Dierentiating(1.5)givesNow,multiplying(1.2)bytheintegratingfactorgives  thatis,ThishastheeectofamelioratingthestilinearpartofthePDE,andwecanuseatime-steppingmethodofourchoice(forexample,afourth-orderRungeKuttafor-mula)toadvancethetransformedequation.Inpractice,onedoesntusetheequationasitiswrittenin(1.8),butratherreplacesactualtime,,withthetimestep,andincrementallyupdatestheformulafromonetimesteptothenext.Thisgreatlyimprovesthestability.Inboththesplitstepmethodandtheintegratingfactormethod,weuseafourth-orderRungeKuttamethodforthetime-stepping.Thefourth-orderRungeKuttaalgorithmthatweusedtoperformthetimeintegrationforthismethodwasc,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-orderAdamsBashforthformula,AM6denotesthesixth-orderAdamsMoultonformula,andAM2denotesamodiedsecond-orderAdamsMoultonformulaspeciedby 23 2Lun1 isthetimestep. A.-K.KASSAMANDL.N.TREFETHENd.Thisequationisexact,andthevariousorderETDschemescomefromhowoneap-proximatestheintegral.IntheirpaperCoxandMatthewsrstpresentasequenceofrecurrenceformulaethatprovidehigherandhigher-orderapproximationsofamulti-steptype.Theyproposeageneratingformulaistheorderofthescheme.Thecoecientsaregivenbytherecurrence 2gm1+1 3gm2+g0 CoxandMatthewsalsoderiveasetofETDmethodsbasedonRungeKuttatime-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.AllenCahnequationwithconstantDirichletboundaryconditions,conditions,1,1],(3.4)u(x,t=0)=47sin(001andthesimulationrunningto=3.Toimposetheboundaryconditionswedeneandworkwithhomogeneousboundaryconditionsinthevariable;thespatialdiscretizationisbyan80-pointChebyshevspectralmethod(seesection5).Weemphasizethattherstthreeproblems,becauseoftheperiodicboundaryconditions,canbereducedtodiagonalformbyFouriertransformation,whereasthefourthcannotbereducedthiswayandthusforcesustoworkwithmatricesratherthanscalars.OurresultsaresummarizedinFigures14.Therstplotineachgurecomparesaccuracyagainststepsizeandthusshouldbereasonablyindependentofmachineandimplementation.Ultimately,ofcourse,itiscomputertimethatmatters,andthisiswhatisdisplayedinthesecondplotineachgure,basedonourMatlabimplementationsonan800MHzPentium3machine.Otherimplementationsandothermachineswouldgivesomewhatdierentresults.BeforeconsideringthedierencesamongmethodsrevealedinFigures14,letusrsthighlightthemostgeneralpoint:ourcomputationsshowthatitisentirelypracticaltosolvethesedicultnonlinearPDEstohighaccuracybyfourth-ordertime-stepping.Mostsimulationsinthepasthavebeenoflowerorderintime,typicallysecondorder,butwebelievethatformostpurposes,afourth-ordermethodissuperior.TurningnowtothedierencesbetweenmethodsrevealedinFigures14,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.4ResultsfortheAllenCahnequation.Thisproblemisnondiagonalandmorechallengingthantheothers.OurMatlabcodeislistedinsectionproblemwithschemesofhigherthansecondorder.Forsecond-ordercalculations,splitstepschemesarecertainlycompetitive.methoddoeswellinallofthediagonalcases.Itisfast,accurate,andverystable.ThisisremarkablewhenwecompareitsperformancewiththatoftheIMEXschemesfromwhichitwasderived.Theonlyproblemwiththisschemeisthedicultyingeneralizingittonondiagonalcases.Finally,theintegratingfactorschemeperformswellfortheBurgersequation.ItdoesntdowellfortheKSequation,though,comingoworstofall,andwecouldntgetittoworkfortheKdVequationortheAllenCahnequationwiththespatialdiscretizationthatweused.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:AllenCahn.TheAllenCahnequationisan-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.9MatlabcodetosolvetheAllenCahnequationandproduceFigure.Again,thiscodetakeslessthansecondtorunonan800MHzPentiummachine. A.-K.KASSAMANDL.N.TREFETHENTREFETHENA.IserlesAFirstCourseintheNumericalAnalysisofDierentialEquations,CambridgeUniversityPress,Cambridge,UK,2000.2000.A.IserlesOnthenumericalquadratureofhighlyoscillatingintegrals:FouriertransformsIMAJ.Numer.Anal.,24(2004),pp.365391.365391.J.C.Jiminez,R.Biscay,C.Mora,andL.M.RodriguezDynamicpropertiesofthelocallin-earizationmethodforinitial-valueproblems,Appl.Math.Comput.,126(2002),pp.6381.6381.G.E.Karniadakis,M.Israeli,andS.A.OrszagHighordersplittingmethodsfortheincompressibleNavier-Stokesequations,J.Comput.Phys.,97(1991),pp.414443.414443.C.A.KennedyandM.H.CarpenterAdditiveRungeKuttaschemesforconvection-diusion-reactionequations,Appl.Numer.Math.,44(2003),pp.139181.139181.J.KimandP.MoinApplicationsofafractionalstepmethodtoincompressibleNavier-Stokesequations,J.Comput.Phys.,59(1985),pp.308323.308323.O.M.Knio,H.N.Najim,andP.S.WyckoffAsemi-implicitnumericalschemeforreacting,J.Comput.Phys.,154(1999),pp.428467.428467.D.KortewegandG.deVriesOnthechangeofformoflongwavesadvancinginarectan-gularcanal,andonanewtypeoflongstationarywaves,Philos.Mag.Ser.5,39(1895),pp.422433.422433.W.KressandB.GustafssonDeferredcorrectionmethodsforinitialvalueboundaryprob-,J.Sci.Comput.,17(2002),pp.241251.241251.M.KrusemeyerDierentialEquations,MacmillanCollegePublishing,NewYork,1994.1994.Y.KuramotoandT.TsuzukiPersistentpropagationofconcentrationwavesindissipativemediafarfromthermalequilibrium,Prog.Theoret.Phys.,55(1976),pp.356369.356369.Y.Maday,A.T.Patera,andE.M.RnquistAnoperator-integration-factorsplittingmethodfortime-dependentproblems:Applicationtoincompressibleuidow,J.Sci.Com-put.,5(1990),pp.263292.263292.R.I.McLachlanandP.AtelaTheaccuracyofsymplecticintegrators,Nonlinearity,5(1992),pp.541562.541562.R.McLachlanSymplecticintegrationofHamiltonianwaveequations,Numer.Math.,66(1994),pp.465492.465492.W.J.MerryfieldandB.ShizgalPropertiesofcollocationthird-derivativeoperators,J.Comput.Phys.,105(1993),pp.182185.182185.P.A.MilewskiandE.G.TabakApseudospectralprocedureforthesolutionofnonlinearwaveequationswithexamplesfromfree-surfaceows,SIAMJ.Sci.Comput.,21(1999),pp.11021114.11021114.M.L.MinionSemi-implicitspectraldeferredcorrectionmethodsforordinarydierentialequa-,Commun.Math.Sci.,1(2003),pp.471500.471500.D.R.Mott,E.S.Oran,andB.vanLeerAquasi-steadystatesolverforthestiordinarydierentialequationsofreactionkinetics,J.Comput.Phys.,164(2000),pp.407428.407428.B.Nicolaenko,B.Scheurer,andT.TemamSomeglobalpropertiesoftheKuramoto-Sivashinskyequation:Nonlinearstabilityandattractors,Phys.D,16(1985),pp.155183.155183.S.P.NrsettAnA-stablemodicationoftheAdamsBashforthmethods,inConferenceonNumericalSolutionofDierentialEquations(Dundee,1969),LectureNotesinMath.109,Springer-Verlag,Berlin,1969,pp.214219.214219.E.OttChaosinDynamicalSystems,CambridgeUniversityPress,Cambridge,UK,1993.1993.R.D.RuthAcanonicalintegrationtechnique,IEEETrans.NuclearScience,NS-30(1983),pp.26692671.26692671.S.J.RuuthImplicit-explicitmethodsforreaction-diusionproblemsinpatternformation,J.Math.Biol.,34(2)(1995),pp.148176.148176.Y.SaadAnalysisofsomeKrylovsubspaceapproximationstothematrixexponentialoperatorSIAMJ.Numer.Anal.,29(1992),pp.209228.209228.J.M.Sanz-SernaandM.P.CalvoNumericalHamiltonianProblems,ChapmanandHall,London,1994.1994.M.SchatzmanTowardnon-commutativenumericalanalysis:HighorderintegrationintimeJ.Sci.Comput.,17(2002),pp.99116.99116.L.M.SmithandF.WaleffeTransferofenergytotwo-dimensionallargescalesinforced,rotatingthree-dimensionalturbulence,Phys.Fluids,11(1999),pp.16081622.16081622.L.M.SmithandF.WaleffeGenerationofslowlargescalesinforcedrotatingstratiedturbulence,J.FluidMech.,451(2002),pp.145169.145169.G.StrangOntheconstructionandcomparisonofdierenceschemes,SIAMJ.Numer.Anal.,5(1968),pp.506517.506517.E.TadmorTheexponentialaccuracyofFourierandChebyshevdierencingmethods,SIAMJ.Numer.Anal.,23(1986),pp.110.