FOURTHORDER TIMESTEPPING FOR STIFF PDEs ALYKHAN KASSAM AND LLOYD N

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-orderRunge�Kuttafor-mula)toadvancethetransformedequation.Inpractice,onedoesn�tusetheequationasitiswrittenin(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,andAM2denotesamodi�edsecond-orderAdams�Moultonformulaspeci�edby 23 2Lun1 isthetimestep. A.-K.KASSAMANDL.N.TREFETHENd
.Thisequationisexact,andthevariousorderETDschemescomefromhowoneap-proximatestheintegral.IntheirpaperCoxandMatthews�rstpresentasequenceofrecurrenceformulaethatprovidehigherandhigher-orderapproximationsofamulti-steptype.Theyproposeageneratingformulaistheorderofthescheme.Thecoe
cientsaregivenbytherecurrence 2gm�1+1 3gm�2+g0 CoxandMatthewsalsoderiveasetofETDmethodsbasedonRunge�Kuttatime-stepping,whichtheycallETDRKschemes.Inthisreportweconsideronlythefourth-orderschemeofthistype,knownasETDRK4.AccordingtoCoxandMatthews,thederivationofthisschemeisnotatallobviousandrequiresasymbolicmanipulationsystem.TheCoxandMatthewsETDRK4formulaeare: an=eLun+L�1(eL�I)N(un), bn=eLun+L�1(eL�I)N(an+h/2), cn=eLan+L�1(eL�I)(2N(bn+h/2)�N(un)), uneLhun+h�2L�3{[�4�Lh+eLh(4�3Lhh)2)]N(un) +2[2+2)+ 4�3Lh�(Lh)2+eLh(4�Lh)]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)andfrom�vetermsofaTaylorexpansion.Forsmall,thedirectformulaisnogoodbecauseofcancellation,buttheTaylorpolynomialisexcellent.Forlarge,thedirectformulais�ne,butthe A.-K.KASSAMANDL.N.TREFETHENThesolutionwehavefoundistoevaluate)viaanintegraloveracontourinthecomplexplanethatenclosesandiswellseparatedfrom0: 2
f(t) becomesamatrixinsteadofascalar,thesameapproachworks,withtheterm1)becomingtheresolventmatrix( HerecanbeanycontourthatenclosestheeigenvaluesofContourintegralsofanalyticfunctions(scalarormatrix)inthecomplexplaneareeasytoevaluatebymeansofthetrapezoidrule,whichconvergesexponentially[15,16,24,60].Inpracticewetaketobeacircleandusually�ndthat32or64equallyspacedpointsaresu
cient.Whenisreal,wecanexploitthesymmetryandevaluateonlyinequallyspacedpointsontheupperhalfofacirclecenteredontherealaxis,thentaketherealpartoftheresult.ThescalarsoreigenvaluesofthatariseinadiscretizedPDEtypicallylieinornearthelefthalfofthecomplexplaneandmaycoverawiderange,whichgrowswiththespatialdiscretizationparameter.Fordiusiveproblemstheyareclosetothenegativerealaxis(e.g.,theKSequation),andfordispersiveproblemstheyareclosetotheimaginaryaxis(KdV).Suitablecontoursmayaccordinglyvaryfromproblemtoproblem.Ourexperienceshowsthatmanydierentchoicesworkwell,solongasoneiscarefultoensurethattheeigenvaluesareindeedenclosedby.Forsomediusiveproblems,itmightbeadvantageoustouseaparaboliccontourextendingto,takingadvantageofexponentialdecaydeepinthelefthalf-plane,butwehavenotusedthisapproachfortheproblemstreatedinthispaper.Fordiagonalproblems,wehavetheadditional�exibilityofbeingabletochooseacontourthatdependson,suchasacirclecenteredat.Inthisspecialcase,thecontourintegralreducessimplytoameanof)over,whichweapproximatetofullaccuracybyameanoverequallyspacedpointsalong(oragainjusttheupperhalfof,followedbytakingtherealpart).Forthedetailsofexactlyhowthiscontourintegralapproachcanbeimplemented,seetheMatlabcodeslistedinsections4and5.Thereisconsiderable�exibilityaboutthisprocedure,andwedonotclaimthatourparticularimplementationsareoptimal,merelythattheyworkandareeasytoprogram.Fornondiagonalproblems,quiteabitofcomputationisinvolved�say,32matrixinverses�butasthisisdonejustoncebeforethetime-steppingbegins(assumingthatthetimestepsareofa�xedsize),theimpactonthetotalcomputingtimeissmall.Itwouldbegreaterforsomeproblemsinmultiplespacedimensions,butinsomecasesonecouldamelioratetheproblemwithapreliminarySchurfactorizationtobringthematrixtotriangularform.ContourintegralsandTaylorseriesarenottheonlysolutionsthathavebeenproposedforthisproblem.BothBeylkin[5]andFriesneretal.[21],forexample,useamethodthatisbasedonscalingandsquaring.Thatmethodisalsoeective,butthecontourintegralmethodappealstousbecauseofitsgreatergeneralityfordealingwitharbitraryfunctions.Todemonstratetheeectivenessofourstabilizationmethod,Table2.2considersthecomputationofthecoe
cient)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.Toimposetheboundaryconditionswede�neandworkwithhomogeneousboundaryconditionsinthevariable;thespatialdiscretizationisbyan80-pointChebyshevspectralmethod(seesection5).Weemphasizethatthe�rstthreeproblems,becauseoftheperiodicboundaryconditions,canbereducedtodiagonalformbyFouriertransformation,whereasthefourthcannotbereducedthiswayandthusforcesustoworkwithmatricesratherthanscalars.OurresultsaresummarizedinFigures1�4.The�rstplotineach�gurecomparesaccuracyagainststepsizeandthusshouldbereasonablyindependentofmachineandimplementation.Ultimately,ofcourse,itiscomputertimethatmatters,andthisiswhatisdisplayedinthesecondplotineach�gure,basedonourMatlabimplementationsonan800MHzPentium3machine.Otherimplementationsandothermachineswouldgivesomewhatdierentresults.BeforeconsideringthedierencesamongmethodsrevealedinFigures1�4,letus�rsthighlightthemostgeneralpoint:ourcomputationsshowthatitisentirelypracticaltosolvethesedi
cultnonlinearPDEstohighaccuracybyfourth-ordertime-stepping.Mostsimulationsinthepasthavebeenoflowerorderintime,typicallysecondorder,butwebelievethatformostpurposes,afourth-ordermethodissuperior.TurningnowtothedierencesbetweenmethodsrevealedinFigures1�4,the�rstthingwenoteisthatthedierencesareveryconsiderable.Themethodsdierine
ciencybyfactorsasgreatas10orhigher.Wewerenotabletomakeeverymethodworkineverycase.Ifamethoddoesnotappearonagraphitmeansthatitdidnotsucceed,seeminglyforreasonsofnonlinearinstability(whichperhapsmighthavebeentackledbydealiasinginourspectraldiscretizations).Akeyfeaturetolookforinthe�rstplotineach�gureistherelativepositioningofthedierentmethods.Schemesthatarefurthertotherightforagivenaccuracytakefewerstepstoachievethataccuracy.Itispossiblethateachstepismorecostly,however,sojustbecauseaschemeachievesagoodaccuracyinfewstepsdoesnotmeanthatitisthemoste
cient.Thesecondplotinthe�guresgivesinsightintothesecomputationtimes.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.Theonlyproblemwiththisschemeisthedi
cultyingeneralizingittonondiagonalcases.Finally,theintegratingfactorschemeperformswellfortheBurgersequation.Itdoesn�tdowellfortheKSequation,though,comingoworstofall,andwecouldn�tgetittoworkfortheKdVequationortheAllen�Cahnequationwiththespatialdiscretizationthatweused.Thisisalsoalittlesurprising,consideringhowwidelyusedthisschemeis. A.-K.KASSAMANDL.N.TREFETHEN Relative timestepRelative error at t = 3 radius = 1e3 Fig.5LossofstabilityintheETDschemeappliedtotheBurgersequationastheradiusofthecontourshrinkstozero.ThecontourofzeroradiuscorrespondstoanETDcalculationdirectlyfromthede�ningformula.Weusetheinitialcondition0)=cos(16)(1+sin(Astheequationisperiodic,wediscretizethespatialpartusingaFourierspectralmethod.TransformingtoFourierspacegives or,inthestandardformof(1.2),u,t denotesthediscreteFouriertransform.WesolvetheproblementirelyinFourierspaceanduseETDRK4time-steppingtosolveto=150.Figure6showstheresult,whichtooklessthan1secondofcomputertimeonourworkstation.Despitetheextraordinarysensitivityofthesolutionatlatertimestoperturbationsintheinitialdata(suchperturbationsareampli�edbyasmuchas10upto=150),wearecon�dentthatthisimageiscorrecttoplottingaccuracy.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.IserlesAFirstCourseintheNumericalAnalysisofDi
erentialEquations,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-di
usion-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.KrusemeyerDi
erentialEquations,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:Applicationtoincompressible�uid�ow,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-surface�ows,SIAMJ.Sci.Comput.,21(1999),pp.1102�1114.1102�1114.M.L.MinionSemi-implicitspectraldeferredcorrectionmethodsforordinarydi
erentialequa-,Commun.Math.Sci.,1(2003),pp.471�500.471�500.D.R.Mott,E.S.Oran,andB.vanLeerAquasi-steadystatesolverforthesti
ordinarydi
erentialequationsofreactionkinetics,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-stablemodi�cationoftheAdams�Bashforthmethods,inConferenceonNumericalSolutionofDi
erentialEquations(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-di
usionproblemsinpatternformation,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.WaleffeGenerationofslowlargescalesinforcedrotatingstrati�edturbulence,J.FluidMech.,451(2002),pp.145�169.145�169.G.StrangOntheconstructionandcomparisonofdi
erenceschemes,SIAMJ.Numer.Anal.,5(1968),pp.506�517.506�517.E.TadmorTheexponentialaccuracyofFourierandChebyshevdi
erencingmethods,SIAMJ.Numer.Anal.,23(1986),pp.1�10.