METHODOFLINESTRANSPOSE:HIGH-ORDERSCHEMESFORPARABOLICPROBLEMSByHanaChoADISSERTATIONSubmittedtoMichiganStateUniversityinpartialentoftherequirementsforthedegreeofAppliedMathematics-DoctorofPhilosophy2016ABSTRACTMETHODOFLINESTRANSPOSE:HIGH-ORDERSCHEMESFORPARABOLICPROBLEMSByHanaChoInthedissertation,wemainlyconsiderdevelopingtnumericalschemesforAllen-CahnandCahn-Hilliardequations,whicharetheoriginoftheldequations.Inthepartofthedisseration,wepresentanewsolverfornonlinearsecond-orderparabolicproblemsthatisL-stableandachieveshighorderaccuracyinspaceandtime.Thesolverisbuiltbyconstructingasingle-dimensionalheatequationsolverthatusesfastO(N)convolution.ThisfundamentalsolverisbasedontheuseoftheGreen'sfunc-tiontoinvertamoHelmholtzequation.Higherordersofaccuracyintimearethenconstructedthroughanoveltechniqueknownassuccessiveconvolution,whichfacilitateourproofsofstabilityandconvergence,andpermitustoconstructschemesthathaveprovabledecay.Themulti-dimensionalsolverisbuiltbyrepeatedapplicationofdimensionallysplitindependentfundamentalsolvers.Wealsosolvenonlinearparabolicproblemsbyusingtheintegratingfactormethod,whereweapplythebasicschemetoinvertlinearterms(thatlooklikeaheatequation),andmakeuseofHerminterpolantstointegratetheremainingnonlinearterms.Oursolverisappliedtoseverallinearandnonlinearequationsincludingheat,Allen-Cahn,andtheFitzhugh-Nagumosystemofequationsinoneandtwodimensions.Inthesecondpartofthedisseration,weextendourMethodOfLinesTranspose(MOLT)schemetoCahnHilliard(CH)andvectorCahnHilliard(VCH)equations.OurstepistoestablishthegradientstabilityforCHmodels.Thisprocedureisjustasimplechangeofvariableswhereoneofthepointsissubtractedfromtheoriginalvariable.Weprovethatinthesemi-analyticsetting,usingBackwardsEulertimestepping.Afterdiscretizingintime,weproceedtoourspatialsolverforinvertingthelinearpartofthesemi-analyticoperatorontothenon-linearparttoconstructantdpointmethod.ThisisdonebyfactoringthefourthintothemoHelmholtzoperators,whichabove.Byincludingthesplittingerrorintotherighthandsideoftheointmethod,wearriveatanon-splitscheme.WealsocombinetheMOLTformulationwithexistingtimesteppingforhighordertimesteppingmethods,andnumericallydemonstratedthegradientstablepropertyin1Dand2Dinallsimulationsrun.Timeadaptivemethodsareshowntobemorentthanusingthesamemethodwithlargetimesteps.ACKNOWLEDGMENTSFirstandforemost,Iwouldliketoexpressmydeepestgratitudetomyadvisor,professorAndrewChristlieb,forhissupportandencouragementduringeyears.Andrewhaspro-videdmewithtremendousknowledge,greatinsight,patientguidanceandpleasantresearchenvironment.Itwassuchanamazingopportunitytodiscusstoourgroupmemberseveryweekandtobuildafriendship.IespeciallythankprofessorMattCausley,whogavemeandourgroupvaluableguidanceaboutMOLTscheme.Hewasalwayswillingtospendhistimetoprovidemeacademicknowledgeandpersonalmentoring.IthanktoprofessorKeithPromislow,forhisenlighteningadvicesthroughhisclassesanddiscussions.Hisstronganalysisofmodelshasbeenvaluable.Ialsothanktherestofcommitteemembers,professorYingdaCheng,forteachingmegreatknowledgeinnumericalanalysis,professorShankerBalasubramaniamandprofessorJosePerea,forvaluablefeedbackandcomments.IalsothanktomyKoreanadvisor,professorJunghoYoon,forhiscontinuoussupportacademicallyandpersonallyfortenyears.IthanktoallmydearfriendsatMSU,EWHA,Youngran,JiguchonandLansingchurch,fortheirhelpandprayers.Ihavenodoubtwewillalwaysbefriendswhereverweare.Iameternallygratefultomyfamily:myparents,myparents-in-law,mysisterEunhye(congratsyourwedding!),mysister-in-law'sfamily,andmygrandparents,fortheirlove,unconditionalsupport,andprayers.Mydeepestgratitudeshouldgotomyhusband,SijungOh,forhisendlesslove.Iwouldn'thavemadeitanyofthisworkwithouthispatience,love,andsupport.Heismybestfriend,mybestteacher,andmyBroomtreeofmylife.Lastbutnotleast,IambeyondwordsinmygratitudebutforeverthanktomyJesus,myLord.TheLordismyshepherd,Ilacknothing.AllglorybelongstomySavior.Amen.ivTABLEOFCONTENTSLISTOFTABLES....................................viiLISTOFFIGURES...................................ixChapter1Introduction...............................11.1method................................11.2Reviewofpreviousworks.............................31.2.1MethodofLines(MOL).........................31.2.1.1Eyre'soperatorsplitting....................41.2.1.2Thestabilizedsemi-implicitSDCmethod..........41.2.1.3Fullyimplicitschemewithconjugategradientiteration...51.2.1.4Exponentialtimeencing(ETD)method........61.2.2MethodofLinesTranspose(MOLT)..................71.2.2.1Fastmultipolemethod(FMM)................81.2.2.2Fouriercontinuation-alternatingdirection(FC-AD)method81.2.2.3Dimensionallysplitting-fastconvolutionmethod.......91.3Outlineofthedissertation............................10Chapter2MOLTforParabolicEquations....................122.1Motivation.....................................122.2Firstorderschemefor1DHeatequation....................132.2.1InversionofthemoHelmholtzoperator..............152.2.2Particularsolution............................182.2.2.1Characteristicdecompositionforparticularsolution.....192.2.2.2Aspatialsolverforparticularsolution............202.2.3Homogeneoussolution..........................212.2.3.1Dirichletboundaryconditions.................222.2.3.2Neumannboundaryconditions................222.2.3.3Periodicboundaryconditions.................232.3Higherorderschemefor1DHeatequation...................232.3.1Convergence................................262.3.2Stability..................................282.3.3Numericalresult:tstudyof1DHeatequation.......332.4MOLTforMulti-DHeatequation........................332.4.1Stability..................................342.4.2Numericalresult:tstudyof2DHeatequation.......352.5MOLTforNonlinearequations.........................362.5.1Allen-CahnEquation...........................382.5.1.1Numericalresult:1DAllen-Cahntest............392.5.1.2Numericalresult:2DAllen-Cahntest............40v2.5.2FitzHugh-NagumoSystem........................422.5.2.1Numericalresult:2DFitzHugh-Nagumotest.........45Chapter3MOLTforCahn-HilliardEquation.................473.1Models.......................................473.1.1transformedmodel............................493.2Firstorderschemefor1DCHequation.....................503.2.1Energystability..............................503.2.2Twononlineariterativeschemes.....................533.2.3Fullydiscretesolution..........................563.3Higherorderschemefor1DCHequation....................593.3.1MOLTwithBackwardFormula(BDF)...........593.3.2MOLTwithSinglyDiagonalImplicitRungeKutta(SDIRK).....603.3.3MOLTwithSpectralDeferredCorrection(SDC)............623.3.4Numericaltest:tstudiesof1DCHsolutions........653.4MOLTforMulti-DCHequation.........................703.4.1Firstorderschemefor2DCHequation.................703.4.2Higherorderschemefor2DCHequation................733.4.3Numericaltest:tstudiesof2DCHsolutions........753.5MOLTforvectorCH(VCH)equation......................763.5.1Twononlineariterativeschemes.....................783.5.2Numericaltest:tstudiesof2DVCHsolutions........793.6Timeadaptivestrategy..............................803.7Numericaltests..................................823.7.11DCahn-HilliardModel.........................823.7.22DCahn-HilliardModel.........................863.7.32DvectorCahn-HilliardModel.....................893.7.42DsixthorderModel...........................92Chapter4ConclusionandFuturework.....................96BIBLIOGRAPHY....................................99viLISTOFTABLESTable2.1:Valuesof2chosenforordersP=1;2;:::6.Thecolumnarethoseusedinourschemes,anduniquelyguaranteedecayandA(0)-stability.Forcomparison,wealsodisplaythevaluesinN˝rsett[49]whichgiveoptimalorderP+1,attheexpenseofdecay..32Table2.2:tstudiesfor1DHeatequation(2.2)((x;t)2(0;2ˇ)[0;4])withPeriodicboundaryconditions.Pindicatestheorderofresolventexpansionin(2.28)............................33Table2.3:tstudyfora2DHeatequation((x;y)2(0;2ˇ)(0;2ˇ))withPeriodicboundaryconditions,(T=1)..............35Table2.4:tstudyforthe1DAllen-Cahn(AC)equationwithaninitialuAC(x;0)in(2.49)andwiththehomogeneousNeumannboundaryconditions.................................39Table2.5:tstudiesfor2DAllen-Cahnequataionwiththeinitial(2.50)andwithhomogeneousNeumannboundaryconditions........43Table3.1:tstudiesofsecond-ordermethodsfor1DCHequationwithperiodicBC................................66Table3.2:tstudiesofthird-ordermethodsforthe1DCHequationwithperiodicBC.............................69Table3.3:tstudyofSDC3(P=3)forthe1DCHequationwithperiodicBC................................69Table3.4:tstudiesofsecond-ordermethodsforthe2DCHequationwithperiodicBC.............................77Table3.5:Ripeningtimeof1DCHequationwithsmalltimestepsizet=0:01)andwiththeperiodicBC..................84Table3.6:Ripeningtimeof1DCHequationwithlargetimestepst=0:05;1;10)andwiththeperiodicBC........................84Table3.7:Ripeningtimeof1DCHequationwithadaptivetimestepsizeandwiththeperiodicBC...........................84viiTable3.8:Numericalripeningtimeofthe2DCHequationwithtimestep-pingmethodst=0:01,Ntol=107)andwiththeperiodicBCs..87Table3.9:Numericalripeningtimeofthe2DCHequation(theperiodicBCs):Comparisonbetweenthetimesteppingmethodst=0:1,Ntol=107)andtheadaptivetimesteppingmethod(BE-BDF2)..89Table4.1:Summarizationofnumericalresultsofpresentedscheme,dimensionalsplitMOLTscheme............................97viiiLISTOFFIGURESFigure2.1:Fouriermodesof(a)theexpnentialoperator,(b)thetruncatedTaylorexpansion(P=6),and(c)theTruncatedLaguerreexpansion(P=6).Thecommonparameters=1andt=0:1ischosen......26Figure2.2:Maximumfactorsj˚(iy)jforthefewordersP,with(a)maximalorder,or(b)decay.Whenmaximizingorder,the3schemesexhibitA-stability,whereasensuringdecayleadstoL-stableschemes.ForP>3,bothschemesbecomeA()-stable.31Figure2.3:Thenumericalsolutionswithtwottimestepsizes,t=0:05andt=0:0063(commonparameters:x=26and=0:03p2).Solid(red)line:theexactin(2.49)atT=8..........40Figure2.4:Timeevolutionofthe2DACnumericalsolutionwithinitialcondi-tion(2.50),uptotimeT=0:02562=10:24.(parameters:t=0:0256;x=y=0:0039;=0:05)..................42Figure2.5:Radiiofthenumericalinterfacialcircleasafunctionoftime(0t0:02562)withvarioustimestepsizes.(a)=0:05(b)=0:01(commonparameter:x=28).Redline:theexactradiusRin(2.52)...................................43Figure2.6:Contourplot:Temporalevolutionoftheconcentrationofactivatoruusingthesecond-orderschemein(2.57)withparameters:t=0:01;x=y=0:2;M=4(fourthorderinspace)..........46Figure3.1:Bifurcationoftheparameter12inthemoHelmholtzoperatorL.....................................55Figure3.2:Energydescentandnonlineariterationcounthistoryofallsecond-andthird-ordermethodswithtimestept=0:025andpa-rameters(3.26)..............................67Figure3.3:Energydescentandnonlineariterationcounthistoryofallsecond-order2Dmethodswithtimestept=0:05andwithparame-ters(3.38).................................77ixFigure3.4:Energydescentandnonlineariterationcounthistoryofallsecond-orderVCH2Dmethodswithtimestept=0:000125andwithparameters(3.45).............................79Figure3.5:Temporalevolutionofthe1DCHsolutionfromtheadaptivetimestep(BE-BDF2)scheme.........................85Figure3.6:Thetimestep,iterationcount,andenergyhistoryofouradaptivetimestep(BE-BDF2)schemefor1DCHequation...........86Figure3.7:Thediscreteenergyoftimesteppingschemes((a)BE,(b)BDF2,and(c)BDF3)for2DCHequation.Thecommonparametersareused:t=5,andx=y=0:05..................87Figure3.8:Iterationcountof2DCHsolutionobtainedbythetimestepp-pingmethod(a)t=0:01),and(b)t=0:1).BlacklineisobtainedbyBackwardEuler(BE)scheme;blueisBDF2;andtheredoneisBDF3atbothplots........................88Figure3.9:Thetimestep,iterationcount,andenergyhistoryofouradaptivetimestep(BE-BDF2)schemefor2DCHequation...........90Figure3.10:Temporalevolutionofthe2DCHsolutionwiththeinitial(3.37)...90Figure3.11:Thetimestep,iterationcount,andenergyhistoryofVCHsolution.93Figure3.12:Temporalevolutionofthe2DCHvectorsolutionwiththeinitial(3.37)foru1and(3.44)foru2.Contoursofcos(argu1+iu2)areplotted.Twophasesu=z2andu=z3havesamecosinevalues,12(lightblueintheplots)butareseparatedbydarkbluelines...93Figure3.13:Theiterationcountandenergyhistoryofthesolutionofsixth-orderequation(3.47)........................95Figure3.14:Temporalevolutionofthe2Dsixthordermodel'ssolutionwiththeinitial(3.37)................................95xChapter1Introduction1.1methodCertainmaterialsthatweuseindailylife,suchasmetals,alloys,ceramicsorpolymers,havepropertiesthatdependontheirmicrostructure.Microstructurecanbeingeneral,ascompositionalandstructuralinhomogeneitiesthatariseduringmaterialprocessing[15].Sp,microstructuresincludegrainarrays,phasedistributions,precipitatedispersions,anddislocationnetworks,eachofwhichproducesdtphysicalphenomena.Sincethemicrostructureofamaterialisstronglyrelatedtoit'sphysicalproperties,composition,orperformance,studyingandinvestigatingitsevolutioniscrucialinmaterialscience.Thisinvestigation,infact,hasledtopracticaluseinmanyareas,suchasmechanics[4],thermodynamics[7],andmorphology[21]andsoforth.Theprocessofmicrostructureevolution,however,isoftenverycomplexduetoenergyinteractions.Sp,thismicrostructureprocessisdrivenbythedecreaseoftotalfreeen-ergy,whichincludeschemicalenergy,interfacialenergy,elasticstrainenergy,electrostaticen-ergyandmagneticenergy[15].Inthiscontext,themodel,whichallowsstudyingthetotalfreeenergy,hasbecomeaconsolidatedtoolforsimulatingmicrostructureevolutionfordecades.Themodelhastwomajorfeatures,interfaceandorderparameter,thusitisworthtoreviewoftheliteratureforinterfaceandorderparameterbr.Aninterfaceisasasurfacewhichformingacommonboundarybetweentwo1tphasesofmatter.Ahistoryofnatureoftheinterfacecanbefoundin[4]andhere,wesummarizetheessenceoftheliterature.Intheearly1800s,theinterfacewasrepresentedbyazerothicknesssurface(sharpinterface)anditwasassumedthatphysicalquantities(e.g.density)werediscontinuousacrosstheinterface.Poisson,MaxwellandGibbsrealizedthattheinterfacerepresentedarapidbutsmoothtransitionofphysicalquantities[30].Inthelate1800s,LordRayleighandVanderWaalsintroducedtheconceptofthenon-zerothicknessinterface(interface),whichcircumventscertaininexplicit-trackinginthesharp-interfacemodel.Theseoriginalconceptsoftheinterfacehavebeenanddevelopedforcenturies.In1935,LandauandLifshitzintroducedaphasetransitioninaferromagneticcrystal[43]suchthatthereisanintermediateregionwherethemagneticmomentschangegraduallyfromonedirectiontotheopposite.Thisistheoriginoftheorderparameter,spaceandtimedependentfuctionwhichiszeroatthedisorderedphaseandnonzeroattheorderedphase.Later,LandauandGinzburgextendedtheideatotheequilibriumofasuperconductor;thefreeenergyofasuperconductorneartransitioncanbedescribedintermsofthespatialvariationofanorderparameter.Landau'sworkhasbeencrucialforthedevelopmentofmoderntheoryofphasetransitionandorderparameterbecameafoundationofmodels.Basedonthisliterature,theoriginoftheequationliesintheworksofCahnandHilliard[8]ontheinterfacialfreeenergyofnon-uniformsystems,andAllenandCahn[3]ontheantiphaseboundarymotioninbinaryalloy.Withoutexplicitlytrackingtheinter-facepositions,bothmodelsareabletopredictcorrespondingcomplexmicrostructure,andthushavebeenextensivelydiscussedinmanyareas,suchassolid-statephasetransformation,graingrowth,andmanyexamplesin[15].2AnimportantgeneralizationofCahn-HilliardmodelisthefunctionalizedCahnHilliard(FCH)model,whichhasbeenproposedbyPromislowet.al.in[29,21,40].TheFCHfreeenergydescribesaphaseseparationinblendsofamphiphilicpolymersandsolvents.More-over,theFCHgradientwsarecomprisedoflong-livednetworkmorphologiesofdistinctco-dimension,andtheauthorsanalyizedtheirgeometricevolution,bifurcationandcompe-tition.AgreatreviewofthedetailsofFCHmodelisin[40].1.2ReviewofpreviousworksInthisSection,wereviewthepreviousworksrelatedtonumericalmethodsforthemodels,especiallyfortheAllen-Cahn(AC)andCahn-Hilliard(CH)equations.Bothequa-tionshavesomenumericalchallenges,suchasthepresenceofthesmallparameterwhichdescribestheinterfacialwidth,thenonlinearityintroducedbythefreeenergy,andvarioustimescalesbetweeneachstageoftemporalevolutions.Therehavebeenmanycontributionstoresolvesuchnumericalchallenges,butwemainlyfocusonfollowingrepresentativemeth-ods:(1)themethodoflines(MOL),and(2)themethodoflinestranspose(MOLT).AreviewofmorediversenumericalapproximationsoftheCHmodelcanbefoundin[56].1.2.1MethodofLines(MOL)TheMethodofLines(MOL)isaclassicaltechnique[52]fortime-dependentpartialtialequations(PDEs).MOLdiscretizesthespatialderivatives,usingvariousmethodssuchasnces,andleavesthetimevariablecontinuous.Thisleadstoasystemofcoupledordinaryerentialequations(ODEs),withthecorrespondinginitialandboundaryconditions.ThebasicideaoftheMOListoreplacethespatial(boundaryvalue)derivatives3inthePDEwithalgebraicapproximations.WewillreviewsomenumericalapproximationsthatarebasedontheMOLliterature.1.2.1.1Eyre'soperatorsplittingIn1998,Eyreproposedaconvexsplittingscheme[26]fortheCHequation(3.1).AfterMOLdiscretization(usingsecond-ordercenteredinspace),hesplitthepotentialfunctionintheCHequation,producingtheconvexandnon-convexfunctionsF(Uj)=Fc(Uj)+Fe(Uj);FcandFearestrictvelyconvex;whereUjapproximatesu(xj;t).Hethentreatedtheconvextermimplicitlyandthenon-convextermexplicitly,f(Unj;Un+1j)=F0c(Un+1j)+F0e(Unj)(Un+1j)3Unj;wheref(u)=F0(u).Thisoperatorsplittingguaranteesunconditionalenergystability,andsolvabilityforanytimestept.ThusmanynumericalschemesformodelsutilizeEyre'ssplitting,suchas[18],whereconvexsplittingiscombinedwithadirectiterativesolver.Fromthepointofviewofsolver,thesplittingmethodsmightbeideal.However,[19,37]indicatethatEyre'ssplittingcanleadtodisproportionatelylargetemporalerrorsduringripening,andcanresultinpoordynamicswhentime-accuratesolutionsarerequired.1.2.1.2Thestabilizedsemi-implicitSDCmethodThestabilizedsemi-implicitspectraldeferredcorrection(stabilizedSISDC)methodwasproposedbyShenet.al.forCHequationandrelatedsystems[46].4Theauthorsusedaspectral-Galerkinmethod(Legendre-Galerkin)inspacewhichcanad-dressnon-periodicboundaryconditions,andappliedthespectraldeferredcorrection(SDC)fortimediscretization,wherethebasicsolver(prediction)isastabilizedsemi-implicitscheme[54]whichistheauthor'spreviouswork.Thereactiontermf(u)wastreatedexplicitly,butanextradissipativetermS2(un+1un)wasintroduced.Toensureunconditionalenergystability,theauthorsalsoconsideredatruncatedpotential~Fsuchthatj~F00(u)jisuniformlyboundedinthemaximumnorm,andchoseSdependingontheupperboundofthej~F00(u)j.Thisapproachguaranteeshighorderunconditionalenergystablemethodsforgeneralgradientwmodels,includingathinmodel,andcomputationalsincethemethoddoesnotrequiresolvinganonlinearsystem.Theadditionalstabilizingtermtheyintroducedinthescheme(prediction),however,leadstotheaccuracylostsothatitmightbeonlyrestrictedtoSDCframeworktogethigherorderofaccuracyintime.Moreover,[18]mentionsthatthissemi-implicitschemeshouldimposeatimesteprestrictiontoensureuniquesolvabilityofCHequation.1.2.1.3FullyimplicitschemewithconjugategradientiterationIn[19],Wettonet.al.presentedanewapproach,afullyimplicittimesteppingschemewithaconjugategradient(CG)spatialsolver,forenergygradientwsfromseveralmodelsincludingAC,CH,higher-orderderivativemodelandvectorvariants.Sp,theyappliedastandardFourierpseudo-spectraldiscretizationinspace,butthisproducedadenseJacobianmatrix,whichin[19]waspreconditionedtlybasedonphysics.TheycouldgethighorderofaccuracyintimeusingtheBackwardformula(BDF)aswellastimeadaptivityforlongtimecoarseningprocess.Theirhighorderschemewasabletoeasilyextendtosixthorderproblemandvectorproblem.5Oneofthemaincontributionsofthisworkisdevelopingseveralbenchmarkproblemsformovinginterfacemodelsin1D,2Dand3D.Inthisdissertation,wewillalsoreproducetheirbenchmarkproblemsofCH,vectorCHandsixthordermodelsinChapter3.However,sincetheirschemeutilizesaFourierspectralmethod,theproblemisrestrictedtotheperiodicboundaryconditions.Moreover,eventhoughtheirnumericalsimulationsshowenergydecayingsolutionforCHmodels,theydonotgiveananalyticalproofforenergystability.1.2.1.4Exponentialtime(ETD)methodTheexponentialtimerencing(ETD)[20],sometimescalledtheexponentialintegrator,isaclassicalnumericalmethodforsystems.TheETDschemeisemployedbytheexactintegrationofthegoverningequationsfollowedbyanapproximationofanintegralinvolvingthenonlinearterms.In[37],aphysics-basedstableschemesforCHandFCHmodelwereproposed.AfteraFourierpseudo-spectraldiscretization,theETDmethodwasemployed.Inparticular,theimplicithighorderRunge-Kutta(IRK)andthebackwardtaylorexpansionswereconsidered.Thisledtoaparallelizablecode,whichwasimplementedonGraphicsProcessingUnits(GPUs).Amaincontributionofthephysics-basedETDschemeisimplementationsofvariousCHandFCHproblems(2D,3D)quicklyandtlyasaparallelcomputing,andtheschemecanpredictvariousgeometriceventsofFCHmodels.However,thephysics-basedETDschemeisalsolimitedtoperiodicboundaryconditionsonrectangulardomains.Anotherissueisthattheschemefromalossofconvergenceorderforlargetimesteps(calledfreezeoutsolution),prohibitingtouseoflargetimesteps6duringcoarsengprocessofmodels.1.2.2MethodofLinesTranspose(MOLT)AsshowninSection1.2.1,severalschemesbasedonMOLdiscretizationhavehadmanycontributionstothenumericalsolutionofmodels.However,sinceMOLschemediscretizesinspacetheboundaryconditionsmustbechosenimmediatelyandanyproofofstabilitydependsonthatdiscretization.TheMOLTscheme,referredtoasRothe'smethod[38],orthehorizontallinemethod,isanalternativetoMOLformulation.TheMOLTstartsbydiscretizingaPDEintimebutcontinuousinspace.Hence,anysta-bilitypropertiesareindependentofthespatialdiscretizationandthePDEistransformedtoacoupledsetofellipticboundaryvalueproblems(BVPs),whichcanbesolvedthroughan-alyticinversion.Thechallenge,then,istoapproximatetheintegralequations.Forinstance,anaiveapproachleadstoadirectsolverformedbyadensematrix-vectorproduct,Un(xj)=ZGfdxˇNXi=1G(xixj)fi;i=1;;N:whereNisnumberofspatialgridpoints,andG()isaGreen'sfunction(kernel).ThenthissummationyieldsO(N2)complexity,whichisimpracticalinlargenumberofN.Duetothelackoftintegralsolvers,theMOLTschemerequiredgreaterstoragethanothermethods,thustheschemedidnotgetagreatdealofattentionbeforetheearly1980s.However,in1987,GreengardandRokhlinsuggestedapotentialtheoryanddevelopedafastalgorithmforintegralequationbasedonthetheory[31],sincethentheMOLTschemehasreceivedmoreattentionrecentlywithcontinuousprogressofthefastalgorithmforintegralequations.7InthisSection,wewillreviewsomerepresentativeworksbasedontheMOLTframework.1.2.2.1Fastmultipolemethod(FMM)Thefastmultipolemethod(FMM)wasintroducedbyGreengardandRokhlintoevalu-atetheCoulombpotentialandforceinparticlesystem[31].TheFMMhassincebeenextendedtoincludeotherpotentialfunctions,suchastheYukawapotential[32],whichisalsocalledaascreenedCoulombpotential.ThebasicideaofFMMisanadaptivequad-treestructureinordertoimposeahierarchyoftsonthecomputationaldomain.Foreachparticle,thenearbyparticlestothepotentialarehandleddirectlybysummation,but(non-neighbor)interactionsarehandledusingmultipoleexpansions[42].Inotherwords,theFMMformsthemulipoleexpansionsforallofnodesinthequad-treestructure,andthenconstructsalocalexpansionforeachnodeinSp,Kropinskiet.al.presentedFMMforthemoHelmholtzequationwhichappearsinthesemi-discretesolutionoftheHeat,theNavier-Stokes,andlinearizedPoisson-Boltzmannequationin[42].BasedonFMM,thepresentedsolutioncomputesasrapidO(N)operations,whereNisthenumberofnodesinthediscretizationoftheboundary,eveninthehighlycomplexdomain.TheFMMisoneofthemosttandrobustmethods,whichcanspeedupthecalculation;however,thisschememightbetformodernparallelprocessingbecauseofglobalcouplingoftheboundaryterms.1.2.2.2Fouriercontinuation-alternatingdirection(FC-AD)methodIn2010,BrunoandLyondevelopedanewmethodforgeneralPDEs,whichisFC-AD(Fourier-ContinuationAlternating-Direction)in[6,48].Thewell-knownAlternatingDirec-8tionImplicit(ADI)approachispioneeredbyDouglas,PeacemanandRachford[22,27].First,theFC(Gram)solveristhemethodofspatialtiation,whichenableshighorder/spectralconvergenceofFourierexpansionsofnon-periodicfunctions.ThehighorderFC-ADalgorithmyieldsunconditionalstabilityforgeneraldomainsatanO(Nlog(N))costpertimestep.TheexactGreen'sfunctionisreplacedwithadimensionallysplitGreensfunction.Asaresult,theboundaryconditionsarenolongergloballycoupled.TheschemealsoenablestparallelimplementationsforgeneraldomainsanddoesnotbothGibbsphe-nomenonandCFLrestrictionsinhyperbolicproblems.However,Thisisattheexpenseofintroducingsplittingerror.1.2.2.3Dimensionallysplitting-fastconvolutionmethodRecently,Christlieb'sgrouphaspresentedthefastconvolutionstrategybasedondimensionalsplittingMOLTschemesin[11,13,10,12].Eachoftheseschemesisabletomakeanalyticallysolvesimpler,one-dimensionalboundaryvalueproblem,andthesubsequentsolutioncanbeconstructedthroughdimensionalsweeps.Thisisalsounconditionallystablenumericalschemewithcomputationalcostandcodingcomplexitycomparabletoexplicitschemes.TheresultingdimensionallysplitMOLTmethodhasbeendemonstratedonarangeofproblemswithcomplexgeometryandboundaryconditions.In[13]successiveconvolutionwasintroduced,leadingtoanA-stablemethodofarbitraryorderintimeforhyperbolicproblems.Moreover,[10,58]includethedevelopmentofanembeddedboundarymethodforNeumannboundaryconditionsonacomplexboundarygeometryandseveralnumericalexamples,suchasMaxwellsequation.Thecontributionsofthepresentauthorto[12]extendthedimensionallysplitMOLTtohighorderinspaceandtimeforlinearandnon-linear9second-orderparabolicproblems,whichwillbereproducedinthisdissertation.Thisworkisongoingtohighorderderivativemodels,suchasCHandFCHequations.Inthesecondpartofthedissertation,weturnourattentiontoabovemodels.Ineachcase,thespatialoperatoriffactored,sothattheboundaryvalueproblemissolvedusingacompositionofone-dimensionalsolves.EachsolverisO(N)andmatrix-free,andthecomplexityiscomparabletoexplicit.1.3OutlineofthedissertationInthedissertation,wemainlyconsiderthetnumericalschemesforAllen-CahnandCahn-Hilliardequations.First,wewouldliketosolvetheequationsinsuchawayastoachievehighorderofaccuracyintimeandspace.Second,wewouldliketomaintaintheenergystabilitywhenthemethodisusedtomodelanymicrostructuralevolution,suchascoarseningprocessinbinaryalloy.InChapter2,weproposeanumericalmethodfortheparabolicequations,includethelinearheatequationandnonlinearequationsoftype.ThisChapterisorganizedasfollows:InSection2.2,wederivethebasicschemefortheone-dimensionalheatequation,whichisL-stableandcanachievehighordersofaccuracyinspaceandtime.InSection2.3,wedescribehowtoobtainanarbitraryorderdiscretizationinasingledimensionwithresolventexpansions.InSection2.4.1,wedescribehowthiscanbeextendedtomultipledimensions,andwepresentresultsforlinearheatinoneandtwodimensionsinSection2.3.3and2.4.2,respectively.InSection2.5,wedescribehowourapproachcanhandlenonlinearsourceterms,andwepresentnumerousnumericalresultsincludingAllen-CahnandtheFitzhugh-Nagumosystemofequations.10InChapter3,weextendthedimensionallysplitMOLTforthe1Dand2DCahn-Hilliardandvectorvariant.InSection3.1,weintroducetheCHmodelandvectorCHmodel,whichareofinterestinthisChapter.InSection3.2,wederiveaorderschemeforCHequationinabasic1Dsettingandtheenergystabilityproofinsemi-discretesettinginSection3.2.1.InSection3.3,wemodifythetraditionaltimesteppingschemetoachievehigherordersofaccuracyintimeandpresentthetemporaltstudiesinSection3.3.4.Inaddition,weextendour1DsolvertomultiplespatialdimensioninSection3.4withthetstudiesin2DinSection3.4.3,aswellasextendtovectormodelinSection3.5.Finally,wedescribeanadaptivetimesteppingstrategyinSection3.6andpresentnumerousnumericalresultsincluding1D,2DandvectormodelinSection3.7.11Chapter2MOLTforParabolicEquations2.1MotivationTheprototypicalparabolictialequationistheheatequation.ut=u;(2.1)subjecttoappropriateinitialandboundarycondition.Heret>0andx2whereˆRnisopen.TheLaplacianistakenwithrespecttothespatialvariablesx=(x1;x2;;xn):u=Pni=1uxixi.Itformsacornerstoneofmathematicsandphysics,anditsunderstandingisessentialformorecomplicatedmathematicalmodels.Fourierintroducedthisequationasameanstodescribetransientheatw.Fickquicklyrecognizeditsimportancetoparticleandchemicalconcentrations.Asaresult,parabolicequationsarenowubiquitousindescribingprocesses,whicharefoundinavastarrayofphysicalproblems,amongwhicharemodelsofchemicalkinetics[3,8],phasemodelsdescribingmorphologyandpatternformationinmultiphaseandsolids[5,9,21],andeventhevolatilityofstocksandbondsinmathematical[55].Numericalsolutionsof(linearandnonlinear)equationshavebeenthesubjectofactiveresearchformanydecades.Asearlyasthe1950'sand60's,itwasrecognizedthatduetotheparabolicscaling,methodoflines(MOL)discretizationsoftheheatequationleadto12numericallysystemsofequations.Inprincipal,thenumericalcanbesubsidedbytakinglargertimesteps,whichareonlystableiffullyimplicitsolversareused.Butinpractice,thememoryofearlycomputerswasextremelylimited,makingfullmatrixinversionsandcostly.Thus,alternatedimensionallyimplicit(ADI)splittingmethods[22,27],whichutilizedimensionalsplittingandtridiagonalsolvers,quicklygainedpopularity.Lateron,memoryconstraintsnolongerdthebottleneckforcomputing,andat-tentionshiftedtowardmethodsthatfocusedonreducinggpointoperations(FLOPs),albeitwithadditionalmemoryconstraints.MostnotableamongtheseareKrylovmethods[35],boundaryintegralmethodsinSections1.2.2.1,andquadraturemethods[47,36].How-ever,withtheadventofGPUprocessors,itappearsthatweareyetagainseeingaparadigmshifttowardsmethodsthatshouldemphasizesmallmemoryfootprints,evenattheexpenseofincurringahigheroperationcount.Thus,ADI-likemethods,whichcantlydecom-poselargerproblemsandlimitoverheadcommunication,warrantfurtherinvestigation,andthesefeaturesarethemotivatingfactorforthiswork.Wenowstarttoproposeanovelnumericalmethod,MOLTinSection2.2forobtainingsolutionstothelinearheatequation(2.1).2.2Firstorderschemefor1DHeatequationWebeginbyformingasemi-discretesolutiontothe1DheatequationusingMOLTscheme.Letu=u(x;t)satisfyut=uxx;(x;t)2(a;b)[0;T];(2.2)13withconstantcot,andappropriateinitialandboundaryconditions.TheMOLTamountstoemployingaschemeforthetimederivative,andcollo-catingtheLaplaciantermattimelevelst=tnandt=tn+1,sothatcollocationhastheformun+1unt=@xxun+un+1un2;>0;whereweintoducedafreeparameter>0.Next,weintroducethetialoperatorcorrespondingtothemoHelmholtzequation,byL=I@xx2;=pt:(2.3)Aftersomealgebra,wethattheschemecanbewrittenasL[un+1(12)un]=2un:(2.4)Wenotethatthereareatleasttworeasonablestrategiesforchoosing:1.Maximizetheorderofaccuracy.Forexample,ifwechoose2=2,thenthediscretizationisthetrapezoidalrule,whichissecondorderaccurateandA-stable.2.Enforcedecay.Forexample,ifwechoose2=1,thenthediscretizationistheBackwardEulerscheme,whichisorderaccurate,L-stable,yetdoesnotmaximizetheorderofaccuracy.Hereandbelow,weoptforthesecondstrategy,asthedecayofnumericalsolutionsoftheheatequationisofparamountimportance.InSection2.3.2,wedevelopthisdiscussioninthecontextofhigherorderschemesthatreliesonacarefulselectionofaswellasrepeated14applicationsofasingleinverseoperator.Uponsolvingequation(2.4)forun+1,wethattheequationfortheupdateisun+1=(12)un+2L1[un];(2.5)thatrequiresinvertingamoHelmholtzoperator.WewillaccomplishthisanalyticallybyusingGreen'sfunction.WeshallL1intheensuingdiscussion.2.2.1InversionofthemoHelmholtzoperatorInthisSection,weseekthefree-spaceGreen'sfunctionG(xjy)asthesolutiontoL[G](x)=(xy);<x;y<1;(2.6)whereisthedeltadistributionandListhetialoperatorin(2.3).WeobtainthesolutionG,whichisstatedinthefollowingtheorem.Theorem2.2.1.ThefreespaceGreen'sfunctionsatisfying(2.6)isgivenbyG(xjy)=2ejxyj;=pt:(2.7)Proof.Weconsiderx6=y,anduponsolvingthehomogeneoustialequation,L[G]=0,G(xjy)=8>>><>>>:Ae(xy)+Be(xy);<x<y;Ce(xy)+De(xy);y<x<1:15ToensurethatGremainsboundedforallxandy,wesetA=D=0.ThenitremainstodetermineBandC.WedemandthatGbecontinuousatx=y,butallowforthederivativeofGtobediscontinuous.Thesizeofthejumpdiscontinuityisobtainedbyintegratingtheequationovertheinterval[yy+]andletting!0,fromwhichwe[Gx]x=y=2:ImposingthatthecontinuityandjumpconditionsyieldsthefollowingsetofequationsCB=0;(C+B)=2:ThesolutionisthengivenbyB=C=2,andtheresult(2.7)follows.WenotethatthemoHelmholtzequationhasanon-oscillatoryGreen'sfunction,andespeciallyitisexponentialdecayingproperty.Now,weconsideradtialequationoftheformL[un+1](x)=un(x);axb;(2.8)whereunisanapproximationofu(x;tn),(n2N).UponmultiplyingbyGandintegratingover[a;b],weafterintegrationbypartsthatun+1(x)=2Zbaejxyjun(y)dy+12u(y)@@yejxyjejxyj@@yu(y)y=by=a:Theboundarytermscanalsobeevaluatedexplicitly,usinglimy!a@@yejxyj=e(xa);limy!b@@yejxyj=e(bx);16whichgives12u(y)@@yejxyjejxyj@@yu(y)y=by=a=Bae(xa)+Bbe(bx);forsomecotsBaandBb.Moreover,theboundarytermssatisfythehomogeneoustialequationL[u]=0.Basedonthisresult,aninversionofthemoHelmholtzoperator(2.3)isthatL1[un]=I@xx21[un]:=2Zbaejxyjun(y)dy+Bae(xa)+Bbe(bx);(2.9)andthecotsBaandBbaredeterminedbyapplyingprescribedboundaryconditionsatx=a;bwhichwedescribeinSection2.2.3.Remark1.Alternatively,hadwefollowedthemethodoflines(MOL)anddiscretized(2.2)inspace,thentheentialoperatorLwouldbereplacedbyanalgebraicoperatorL,andwouldbeinvertednumerically.Remark2.Althoughtheupdate(2.5)(with6=1)isonlyorderaccurate,wedescribeinSection2.3howtoextendourproceduretoarbitraryorderintime.Remark3.Sincedimensionalsplittingisused,allspatialquantitiesarecomputedaccordingtoaone-dimensionalconvolutionintegraloftheform(2.9),whichisperformedonaline-by-linebasis,followingso-called"dimensionalsweeps".SincethediscreteconvolutioniscomputedinO(N)complexity,thefullsolverscaleslinearlyinthenumberofspatialpoints(assumingeachsweepisperformedinparallel).Thus,wethattheinversionoperatorconsistsoftwoparts,theparticularsolutionandthehomogeneoussolution.Theconvolutionoperatoriscomprisedofaparticularsolution,17whichisbytheconvolutionintegralI[un](x):=2Zbaejxyjun(y)dy;n2N;(2.10)andahomogeneoussolutionBae(xa)+Bbe(bx);(2.11)bothofwhichcanbeconstructedinO(N)operationsusingfastconvolution.Furthremore,afullydiscreteschemeisobtainedafteraspatialdiscretizationof(2.9).Thedomain[a;b]ispartitionedintoNsubdomains[xj1;xj],witha=x0<x1<:::xN=b.Wenowdescribeeachoftheseinturn,startingwiththe2.2.2ParticularsolutionTheparticularsolutionissplitintoI[u](x)=IL(x)+IR(x),whereIL(x)=2Zxae(xy)u(y)dy;IR(x)=2Zbxe(yx)u(y)dy:Wenowmakeafewkeyobservationsabouttheparticularsolution,whichwillbeusedextensivelyintheensuingdiscussion.182.2.2.1CharacteristicdecompositionforparticularsolutionEachofthesepartsindependentlysatisfytheorder"initialvalueproblems"(IL)0(y)+IL(y)=2u(y);a<y<x;IL(a)=0;(IR)0(y)IR(y)=2u(y);x<y<b;IR(b)=0;wheretheprimedenotesspatialtiation.Fromtheintegratingfactormethod,theintegralthefollowingidentity,knownasexponentialrecursionIL(xj)=ejIL(xj1)+JL(xj);whereJL(xj)=2Zxjxj1e(xjy)u(y)dy;(2.12)IR(xj)=ej+1IR(xj+1)+JR(xj);whereJR(xj)=2Zxj+1xje(yxj)u(y)dy;(2.13)andj=hj;hj=xjxj1:Thisexpressionisstillexact,andindicatesthatonlythe"local"integralsJLandJRneedstobeapproximated.Bysymmetry,theschemeforJRfollowsfromthatofJL,whichwedescribe.192.2.2.2AspatialsolverforparticularsolutionWeconsideraprojectionofu(y)ontoPM,thespaceofpolynomialsofdegreeM.Alocalapproximationu(xjzhj)ˇpj(z);z2[0;1];isaccuratetoO(hMj),andaquadratureoftheformJL(xj)=j2Z10ejzu(xjhjz)dzˇj2Z10ejzpj(z)dz:(2.14)IfstandardLagrangeinterpolationisused,thenthepolynomialscanbefactorizedaspj(z)=rXk=`pjk(z)uj+k=zTA1juMj;(2.15)wherez=[1;z;:::;zM]T,anduMj=[uj`;:::;uj;:::;uj+r]T,andAjistheVandermondematrixcorrespondingtothepointsxj+k,fork=`:::r,Aj=266666664xjxjlhj0:::xjxjlhjM.........xjxj+rhj0:::xjxj+rhjM3777777752RM+1;M+1Thevaluesof`andraresuchthat`+r=M,andarecenteredaboutjexceptneartheboundaries,whereaone-sidedstencilisrequired.Substitutingthisfactorizationinto(2.14)andintegratingagainstanexponential,wethatJL(xj)ˇJLj:=rXk=`wjkuj+k;20wheretheweightswj=[wj;`;:::wj;r]satisfywTj=˚TjA1j(2.16)andwhere˚jk=j2Z10ejzzkdz=ej2+k!2kj0@1ejk1Xp=0(j)pp!1A(2.17)Iftheweightsarepre-computed,thenthefastconvolutionalgorithmscalesasO(MN)pertimestep,andachievesaO(M)inspace.Ineveryexampleshowninthiswork,wechooseM=2;4orM=6.2.2.3HomogeneoussolutionThehomogeneoussolutionin(2.11)isusedtoenforceboundaryconditions.Weobservethatalldependenceonxintheconvolutionintegral,I(x):=I[un](x),in(2.10)isontheGreen'sfunction,whichisasimpleexponentialfunction.Throughdirecttiation,weobtainIx(a)=I(a);Ix(b)=I(b):(2.18)Variousboundaryconditionsatx=aandx=bcanbeenforcedbysolvingasimple22systemfortheunknownsBaandBb.Remark4.Thecasesofapplyingentboundaryconditionsatx=aandx=barenotconsideredhere,butthedetailsfollowfromananalogousproceduretothatdemonstratedbelow.212.2.3.1DirichletboundaryconditionsLetusbeginwithDirichletboundaryconditions,whichweshallbyun(a)=waandun(b)=wb(ateachdiscretetimestep,t=tn).Thenwa=I(a)+Ba+Bbwb=I(b)+Ba+Bb;where=e(ba).ThissystemofequationscanbesolvedfortheBaandBb,yieldingBa=(waIa)(wbIb)12;Bb=(wbIb)(waIa)12;(2.19)whereIa=I(a)andIb=I(b).2.2.3.2NeumannboundaryconditionsToapplytheNeumannboundaryconditions,weusetheidentities(2.18),anddenoteunx(a)=w0aandunx(b)=w0b.Thisyieldsw0a=I0(a)Ba+Bb;w0b=I0(b)Ba+Bb;andthissystemofequationscanbesolvedfortheBaandBb,Ba=(I0aw0a)+(w0b+I0b)12;Bb=(w0b+I0b)+(I0aw0a)12:(2.20)222.2.3.3PeriodicboundaryconditionsWecanalsoimposePeriodicboundaryconditions,byassumingthatun(a)=un(b);unx(a)=unx(b);8n2N:(2.21)Wenextenforcethisassumptiontoholdonthescheme(2.9),L1[un](a)=L1[un](b)()I(a)+Ba+Bb=I(b)+Ba+Bb;L1x[un](a)=L1x[un](b)()(I(a)Ba+Bb)=(I(b)Ba+Bb);wheretheidentities(2.18)areusedtoL1x.SolvingthislinearsystemyieldsBa=Ib1;Bb=Ia1:(2.22)2.3Higherorderschemefor1DHeatequationInourgroup'srecentwork[13],Causleyet.al.appliedasuccessiveconvolutionapproachtoderivehighorderA-stablesolversforthewaveequation.Weextendthissuccessiveconvolutionideabasedonnewapproach,resolventexpansion,fortheheatequation.Thekeyideaistorecognizethefactthat,inviewofthemoHelmholtzoperator(2.3),thesecondderivativecanbefactoredas@xx2=LI=LIL1:=LD;(2.23)23whereD=IL1;L=(ID)1:(2.24)Substitutionofthesecondexpressioninto(2.23)determinesthesecondderivativecompletelyintermsofthisnewoperator@xx2=(ID)1D=1Xp=0Dp:(2.25)Thisshowsthatsecondorderpartialderivativesofatlysmoothfunctionu(x)canbeapproximatedbytruncatingaresolventexpansionbasedonsuccessivelyapplyingDtou(x),whichisalinearcombinationofsuccessiveconvolutions(2.9).Wealsonotethatforimplementation,eachoperatorisappliedsuccessively,andisbyD(p+1)[u]:=D[Dp[u]];D0[u]:=u:Now,weconsiderasolutionu(x;t)totheheatequation(2.2),thatforsimplicitywetaketobesmooth.WeperformaTaylorexpansiononu(x;t+t),andthenusetheCauchy-Kovalevskayaprocedure[53]toexchangetemporalandspatialderivativestoyieldu(x;t+t)=1Xp=0t@t)pp!u(x;t)=1Xp=0(t@xx)pp!u(x;t)=:et@xxu(x;t):(2.26)Thetermet@xxisaspatialtialoperator,anditcompactlyexpressesthefullTaylorseries.Ourgoalistomakeuseoftheformula(2.25)toconverttheTaylorseriesintoaresolventexpansion[24,25].Thiscanbeperformedterm-by-term,andrequiresrearrangingadoublysum.However,ifweinsteadworkdirectlywiththeoperator24theTaylorseries,thenet@xx=e2@xx2=e2(I)1D:Thelastexponentialformoftheoperatorremindsusageneratingfunctionoforthogo-nalpolynomials,whichisthegeneratingfunctionofthegeneralizedLaguerrepolynomialsL()p(z),1Xp=0L()p(z)tp=1(1t)+1etz1t;z>0;(2.27)andthisbearsastrikingresemblancetoourexpansion.Indeed,ifwetake=1,substitutez=2andt=D,thene2(I)1D=1Xp=0L(1)p(2)Dp=I+1Xp=1L(1)p(2)Dp:(2.28)Therefore,theLaguerreexpansionsviaoursuccessiveconvolutionsofDoperatorandtheTaylorexpansionareexactlysamewithourexponentialformofoperator,e2(I)1D,inthelimitingsense.However,wedoprefertheLaguerreexpansionforthestabilityreason.TocomparetheLaguerreexpansionwiththeTaylorexpansion,weconsiderthetraditionalFourierstabilityanalysis.BeforereplacingtheminaFouriermode,weshouldtruncateeachseriesatsomeorderinordertodevelopanumericalscheme.Thus,weassumethattruncateeachseriesatthesameorderP,whichgivesusTaylor:PXp=0(k2)pp!;Laguerre:PXp=0Lp(2)^Dp;^D=111+k2;wherekisthewavenumber.Then,basedonthepreviousidentities(2.26)and(2.28),25(a)etk2(b)Taylor(P=6)(c)Laguerre(P=6)Figure2.1:Fouriermodesof(a)theexpnentialoperator,(b)thetruncatedTaylorexpansion(P=6),and(c)theTruncatedLaguerreexpansion(P=6).Thecommonparameters=1andt=0:1ischosen.twopartialsumswillconvergetothesamesolutionetk2,thatisexponentialdecayingfunctionintheFouriermode.WespeciallychoseP=6anddraweachFouriermodeoftheexpansionsinFigure2.1.AsshowninFigure2.1b,thetruncatedTaylorexpansiondivergesforlargerwavenumberk,however,thetruncatedLaguerreexpansionatthesameorderPisnicelydecaying,whichasshowninFigure2.1c.ThisisbecausetheFouriermodeofconvolutionoperatortheboundedness,j^Dj1,whateverthewavenumberkis.BecauseofthisboundednessofourDoperator,wedoprefertheLaguerreexpansionforourhigherordertimesteppingscheme.Furthermore,thishasatheoreticalproofabouttherateofconvergenceofpartialsum,whichwillbediscussedinthefollowingSection.2.3.1ConvergenceThisexpansionhasbeenconsideredinthecontextofC0-semigroups[33,1],where@xx2isreplacedwithageneralclosedoperatorAonaHilbertspaceX.Inournotation,werestatepart(ii)ofTheorem5.2in[1],whichisproventherein.26Theorem2.3.1.LettheC0-semigroupT(2)=e2@xx2=1Xp=0L(1)p(2)Dp=1Xp=0L(1)p(2)Lp(LI)p;=1Xp=0L(1)p(2)I@xx2p@xx2p;whereL(1)paregeneralizedLaguerrepolynomials.ThiscanbeapproximatedbyPthpartialsumoftheseries,TP(2)=PXp=0L(1)p(2)Dp:Then,foru2C2P+2,thereexistsforeach2>0anintegerm0suchthatforallintegers2kP+1,withPm0,T(2)uTP(2)uc(2;k)Pk=21@xx2ku;wherec(2;k)isaconstantthatdependsonlyon2andk.Remark5.Thesalientpointofthetheoremisthat,inconsiderationof(c.f.Eqn.(2.3)),theapproximationerrorisoftheformctP+1u(2P+2)(x),whichmatchestheformgivenbyatypicalTaylormethod.Finally,wetruncatetheresolventexpansion(2.28)atorderp=P.Fortheheatequation,thisthenumericalmethodasu(x;t+t)=u(x;t)+PXp=1L(1)p(2)Dp[u](x;t);(2.29)27whichhasatruncationerroroftheform˝:=L(1)P+1(2)DP+1[u](x;t)+OtP+2):(2.30)ForP=1;2;3,theseschemes(evaluatedatt=tn)areun+1=un2D[un];(2.31)un+1=un2D[un]242D2[un];(2.32)un+1=un2D[un]242D2[un]24+66D3[un];(2.33)whereweusedtheRodriguesformulaforgeneralizedLaguerrepolynomialsL(1)0(2)=1;L(1)1(2)=2;L(1)p(2)=zezp!dpdzpezzp1z=2;p>1:Now,thenextnaturalquestionbecomeshowwecanchoosetheparameterforourtimesteppingmethod.Thisparameteractuallyguarenteesthedecayofournumericalscheme.ThiswillbediscussedinthenextSection.2.3.2StabilityThereremainsonecriticalissue:thechoiceofthefreeparameter.In1974,N˝rsettstudiedasimilarsingle-stepmultiderivativemethodfortheheatequation[49]andhetoo,hadafreeparameterinhissolver.WefollowhisleadontheVon-NeumannanalysisbasedonhisMOL28discretization,butinthisworkweoptimizetoobtaindecay,whereasN˝rsettchosetomaximizetheorderofaccuracyofthesolver.Considerthelineartestproblemdydt=;y(0)=1;2C;whoseexactsolutiony(t)y(tn+t)=ezy(tn)1Xp=0zpp!y(tn);z=t2C:Ontheotherhand,ourapproximationsolutionyn+1isobtainedbyPthpartialsumoftheresolventseries(2.28).Applicationof(2.29)tothistestproblemresultsinyn+1=PXp=0L(1)p(2)Dp(z;2)yn;D(z;2)=11z21=zz2:Inotherwords,wecanexpressourapproximationasyn+1=˚P(z)yn=˚P(z)2yn1==˚P(z)n+1;where˚P(z)isthestabilityfunctionofthepresentedPth-ordertimesteppingmethod,˚P(z)=PXp=0L(1)p(2)Dp(z;2);(2.34)whichdependsonthefreeparameter(determinedlater)andthetemporalorderP.Followingstandardions,wesaythatanumericalschemeisA-stable,provided29j˚j1intheleft-halfofthecomplexplanez2C.Likewise,aschemeexhibitsdecayif˚(z)!0asRe(z)!.IfanA-stablemethodalsoexhibitsdecay,itisL-stable.WenotethefollowingrecurrencerelationofgeneralizedLaguerrepolynomials,inordertocalculatethelimitofstabilityfunctionsatnegativenity.ThegeneralizedLaguerrepolynomialssatisfymanyidentities,thefollowingofwhichisthemostpertinent:L(0)p+1(x)L(0)p(x)=L(1)p+1(x)=xp+1ddxL(0)p+1(x):(2.35)Here,L(0)p(x)isthestandardLaguerrepolynomialLp(x).Now,observingthatD!1asRe(z)!,wendthatlimRe(z)˚P(z)=PXp=0L(1)p(2)=L00(2)+PXp=1L0p(2)L0p1(2)=L0P(2);(2.36)wherewehaveusedthetwoexpressionsin(2.35)tointroduceatelescopingsum.Wearenowpreparedtoprovethefollowing.Theorem2.3.2.Letu(x;t)beanapproximatesolutiontotheheatequation(2.2),givenbythesuccessiveconvolutionscheme(2.29).Then,1.If2=x(P)1ischosenasthesmallestrootofL0P+1(x)=(L(0)P+1(x))0,thentheschemeachievesorderP+1,butdoesnotexhibitdecay.2.If2=x(P)1ischosenasthesmallestrootofLP(x)=L(0)P(x),thentheschemeachievesorderP,andexhibitsdecay.3.Followingthestrategy,theschemesareA-stableforP=1;2;3,whereasthesecondstrategyensuresL-stability.Forbothstrategies,A()-stabilityisachievedforP>3,30(a)MaximalOrder(b)DecayFigure2.2:Maximumfactorsj˚(iy)jforthefewordersP,with(a)maximalorder,or(b)decay.Whenmaximizingorder,the3schemesexhibitA-stability,whereasensuringdecayleadstoL-stableschemes.ForP>3,bothschemesbecomeA()-stable.withlargevaluesofˇˇ=2.Proof.Theprooffollowsbyapplyingthemaximummodulusprinciplecoupledwith(2.36).Forpart1,uponexaminingthetruncationerror(2.30),weseethatanadditionalorderofaccuracycanbegainedifwechooseL(1)P+1(2)=2P+1L0P+1(2)=0:However,LP(2)6=0forthischoice,andsodecaydoesnothold.Forpart2,weinsteadenforcedecay,butthenthetruncationerrorisoforderP.Finally,part3isdemonstratedbythemaximumfactors˚alongtheimaginaryaxis,asshownforbothstrategiesinFigure2.2.Inparticular,weobservethatj˚P(iy)j1forP=1;2;3.31DecayMaximalOrderP2LP(2)2LP(2)11.000002.0000-1.000020.585801.2679-0.732030.415800.9358-0.630440.322500.7433-0.576850.263600.6170-0.543660.222800.5277-0.5211Table2.1:Valuesof2chosenforordersP=1;2;:::6.Thecolumnarethoseusedinourschemes,anduniquelyguaranteedecayandA(0)-stability.Forcomparison,wealsodisplaythevaluesinN˝rsett[49]whichgiveoptimalorderP+1,attheexpenseofdecay.Sp,theexplicitformsofstabilityfunctionsare˚1=12zz2;˚2=122+4=2z2+422z+4z22;˚3=132+324166z3+32+34126z2+(346)z6z23:Remark6.In[49],theschemewaschosentomaximizetheorderofaccuracy,implicitlyleadingtoeliminatingtheterminthetruncationerror(2.30),whichisequivalenttothestrategy.However,inthisworkwefollowthesecondstrategy,andchoose2asthesmallestrootofLP(x)toensuredecay.Forcomparisonwerecordthevaluesof2chosenforeachorder1P6,tothoseofN˝rsettinTable2.1.Forallofoursolvers,wechoosetobethelargestpossiblevaluethatstillyieldsprovabledecay.32P=1P=2P=3tL1errororderL1errororderL1errororder0:11:84051041:62551062:42251080:059:21211050:99854:08411071:99283:06201092:98390:0254:60841050:99931:02361071:99643:850110102:99150:01252:30481050:99962:56221081:99824:840210112:99180:006251:15251050:99986:40971091:99906:202110122:9642Table2.2:tstudiesfor1DHeatequation(2.2)((x;t)2(0;2ˇ)[0;4])withPe-riodicboundaryconditions.Pindicatestheorderofresolventexpansionin(2.28).2.3.3Numericalresult:tstudyof1DHeatequationWeillustratetheaccuracyofourmethodforthe1Dheatequationin(2.2).Weconsiderinitialconditionsu(x;0)=sin(x),forx2[0;2ˇ]withperiodicboundaryconditions.WeintegrateuptoatimeofT=4,andset=0:182.Weusethefastconvolutionalgorithmthatisfourthorderaccurateinspace(M=4),andsetthespatialgridsizetobex=2ˇ1024ˇ0:0061.Thisensuresthatthedominanterrorinthesolutionistemporal.WecomputeerrorsbytheL1-norm,andcompareagainsttheexactsolutionu(x;T)=eTu0(x)eTsin(x):(2.37)TheresultofatemporaltstudyforP=1;2and3ispresentedinTable2.2,computedatT=4.Thisexampleservestoillustratethatourschemeachievesthedesignedorderofaccuracyafterappropriatetruncationsoftheresolventexpansions.2.4MOLTforMulti-DHeatequationWeextendthe1Dsolvertomultiplespatialdimensionsthroughtheuseofdimensionalsplitting.Ourkeyobservationisthatwecanusethefactorizationpropertyoftheexponential33toperformtheseriesexpansion.Forinstance,inthreedimensions,wehaveetr2=e2@xx2e2@yy2e2@zz2:(2.38)Now,wereplaceeachtermwiththeidentity(2.28)dimensionbydimension,andthentruncatetheexpansionswhichwillbeintermsoftheunivariateoperatorsL1andDfor=fx;y;zgasdby(2.3),and(2.29)actingonafunctionun(x;y;z).ThissumwiththreeindicesmustthenbetruncatedtoorderP,andafterachangeofindicesweEP=PXp=0PpXq=0P(p+q)Xr=0L(1)p(2)L(1)q(2)L(1)r(2)DpxDqyDrz;(2.39)in3D,withthecorresponding2DoperatorgivenbyEP=PXp=0PpXq=0L(1)p(2)L(1)q(2)DpxDqy:(2.40)Hereweadoptthatsumsaretakenoverallnon-negativeindicesthatsumtoP.2.4.1StabilityTheproofofstabilityforthemulti-dimensionalalgorithmfollowsdirectlyfromthatoftheone-dimensionalcase,withthesameapproachappliedtoeachspatialdimension.Forin-34P=1P=2P=3tL1-errororderL1-errororderL1-errororder0:19:81821058:67171071:29251080:054:91431050:99852:17881071:99281:63541092:98250:0252:45841050:99925:46081081:99632:079110102:97560:01251:22951050:99961:36721081:99792:920410112:8317Table2.3:tstudyfora2DHeatequation((x;y)2(0;2ˇ)(0;2ˇ))withPeriodicboundaryconditions,(T=1).stance,forthe2Dcase,thestabilityfunctioncanbeas˚(z)=PXp=0PpXq=0L(1)p(2)L(1)q(2)DpxDqy;Dx=Dy=11z21;thus,similarwiththe1Dstabilityfunction,limRe(z)˚(z)=L0P(2)2;wherewehaveusedthatDx,Dyapproachto1asRe(z)!andtelescopingsumofLaguerrepolynomial.3Dcaseisalsosimilarwithlowspatialdimensioncase.2.4.2Numericalresult:tstudyof2DHeatequationAsasecondexample,wepresentresultsforthe2Dheatequation.Weconsiderinitialconditionsu(x;y;0)=sin(x)sin(y),for(x;y)2[0;2ˇ][0;2ˇ]withperiodicboundaryconditions.Weuseauniformmeshofsizex=y=2ˇ=512ˇ0:0123.Likewise,theL1-erroriscomputedbycomparingagainsttheexactsolutionu(x;y;T)=e2Tu0(x;y)atthetimeT=1.InTable2.3,wepresentresultsforatemporaltstudyforordersP=1;2;and3.352.5MOLTforNonlinearequationsWenextextendourmethodtononlinearreactiosystemsoftheformut=Dr2u+F(u);(x;t)2(0;T];(2.41)whereu=(u1;u2;;uN),withui=ui(x;t),Disausioncotmatrix,andthereactiontermF:=(f1;f2;;fN)isafunctionofui,(i=1;2;;N).Intheabove,ˆRNisaboundeddomain,andweassumeappropriateinitialvaluesandboundaryconditions.Weshallviewtheasbeingthelinearpartofthetialoperator,andinvertthislinearpartanalytically,usingsuccessiveconvolution.Toderivethescheme,weuseoperatorcalculustorstwrite@tDr2u=F=)eDtr2ut=eDtr2F;(2.42)whereeDtr2isatialoperator.Uponintegrating(2.42)overtheinterval[t;t+t],wearriveattheupdateequationu(t+t)eDtr2u(t)=ZttteD(tt˝)r2F(˝)d˝=Zt0eDt˝)r2F(t+˝)d˝;(2.43)wherewehavemadeuseoftheabbreviatednotation,F(t):=F(u(x;t)).Onthelefthandside,thetermshavebeencollectedbythisentialoperator,andwillbeapproximatedusingthesuccessiveconvolutiontechniquesdevelopedabove.Thereactiontermsontherighthandside(2.43)arefullynonlinear,andwemustconsidernonlinear36stabilitywhenchoosingamethodofdiscretization.Weconsiderapproximatingtheintegralontherighthandside(2.43)withthetrape-zoidalrule.Thisasingle-stepupdateequation,whichwillbesecondorderaccurateu(t+t)eDtr2u(t)=t2eDtr2F(t)+F(t+t):(2.44)Wemayalsoobtainasingle-stepthirdorderscheme,usingmultiderivativeintegration[53].Byreplacingtheintegrand(2.43)withathirdorderinterpolantandper-formingexactintegrationoftheresultingfunction,wearriveatu(t+t)eDtr2u(t)=eDtr2t3F(t)+t6Dtr2F(t)+tdFdt(t)+t3F(t+t);(2.45)wheredFdt(t)=dFdu(t)(Dr2u(t)+F(t)).Theinterpolantthatmatchestheintegrandin(2.43)attimes˝=0,and˝=t,aswellasitsderivativeattime˝=0producesthequadraturerulein(2.45).Uponperusingthethirdorderupdateequation(2.45),wewillneedtousesuccessiveconvolutiontoreplacethetialoperatorexpDtr2,aswellastheLaplacianoperatorr2.Thislatterpointhasbeendetailedin[13],andsowecommentonithere.Usingtheone-dimensionalexpansion(2.25),weobservethatthetwo-dimensionalLaplacianissimilarlygivenbyr22=@xx2@yy2=1Xp=1Dpx+Dpy;andcanbetruncatedattheappropriateaccuracyp=P.Here,thesubscriptsindicatethattheconvolutionisonlyinonespatialdirection,andtheothervariableisheldThus,37Dxisappliedalonghorizontallinesfory-values,andlikewiseforDy.Remark7.Theproposedschemesin(2.44)and(2.45)producenonlinearequationsforu(x;t+t)thatneedtobesolvedateachtimestep.Therefore,anyimplicitsolverwillnecessarilybeproblemdependent.Fortheproblemsexaminedinthiswork,wemakeuseofsimplexed-pointiterativeschemes.WestabilizeouriterativesolversbylinearizingF(u)aboutabackgroundstateFu(u),whichdependsontheproblemunderconsideration.2.5.1Allen-CahnEquationWeexamineingreaterdetailtheapplicationofourschemtotheAllen-Cahn(AC)equation[3],ut=2r2u+f(u);(x;t)2(0;T];(2.46)wherethereactiontermisf(u)=uu3,andˆRdisaboundeddomain,anduhomogeneousNeumannboundaryconditions.Forourpointiteration,welinearizefaboutthestablepointsu=1,whichsatisfyf0(u)=0.Forexample,thesecondorderschemefrom(2.44)becomes(1+t)un+1;k+1=e2tr2un+t2fn+t2fn+1;k+2un+1;k;(2.47)wherenindicatesthetimestepasbefore,andnowkistheiterationindex.Bylaggingthenonlineartermfn+1;k,thepointupdateismadeexplicit.Likewise,thethirdorderschemefrom(2.45)becomes38P=1P=2P=3tL1errororderL1errororderL1errororder0:0252:82161041:38951052:60601060:01251:44191040:96863:61151061:94393:94171072:72490:00637:28741050:98459:21641071:97035:50101082:84110:00313:66321050:99232:32941071:98427:31221092:91130:00161:83651050:99615:86951081:98869:571410102:9335Table2.4:tstudyforthe1DAllen-Cahn(AC)equationwithaninitialuAC(x;0)in(2.49)andwiththehomogeneousNeumannboundaryconditions.1+t3un+1;k+1=e2tr2un+t3fn+t62tr2fn+tfnt+t3fn+1;k+2un+1;k:(2.48)Here,e2tr2isagainunderstoodbyreplacingitwitharesolventexpansion,whichisatruncatedseriesofsuccessiveconvolutionoperators.2.5.1.1Numericalresult:1DAllen-CahntestWedemonstratetheaccuracyofourproposedschemesbysimulatingawell-knowntravelingwavesolution[14,54,44],uAC(x;t)=121tanhxTsst2p2;x2=[0;4];0tT:(2.49)Here,s=3p2=0:09isthespeedofthetravelingwave,andwechoose=0:03p2.Additionally,wechoosethedelaytimeTs:=1:5sT,sothattheexactsolutionuAC(1:5;T)=0:5.ResultsforatimeofT=8areshowninFigure2.3,withtwottimesteps.Thesolutionsagreewellwiththeexactsolution.39Figure2.3:Thenumericalsolutionswithtwottimestepsizes,t=0:05andt=0:0063(commonparameters:x=26and=0:03p2).Solid(red)line:theexactin(2.49)atT=8.InTable2.4,wepresenttheL1-errorinthenumericalsolutionataltimeT=1,usingtheexactsolutionuAC(x;T)(2.49).WeobserveorderaccuracyfromtheBackwardEulermethod,andtheexpectedordersofaccuracyfromthesecond(2.47)andthird(2.48)orderschemes.Toensurethatthetemporalerrorisdominant,wehaveusedthefourthorderaccuratescheme(eq.(2.14)withM=4),withx=29toperformspatialintegrationinthesuccessiveconvolutions.Inprinciple,wecanachievehigherordersaccuracyinspaceandtime.Thelatterwouldrequireusinghigherorderinterpolationtodiscretizethereactiontermin(2.43).2.5.1.2Numericalresult:2DAllen-CahntestWenextsolvetheAllen-Cahnequationintwospatialdimensions.Astandardbenchmarkprobleminvolvesthemotionofacircularinterface[14,54,44],towhichanexactsolutionis40knowninthelimitingcase!0.Theradiallysymmetricinitialconditionsarebyu(x;y;0)=tanh 0:25p(x0:5)2+(y0:5)2p2!;(2.50)whichhasaninterfacialcircle(u(x;y;0)=0)centeredat(0:5;0:5),witharadiusofR0=0:25.Thisinterfacialcircleisunstable,andwillshrinkovertime,asdeterminedbythemeancurvature[3],V=dRdt=1R:(2.51)HereVisthevelocityofthemovinginterface,andRistheradiusoftheinterfacialcircleattimet(i.e.,itistheradiusofthecurvebyu(x;y;t)=0).Byintegrating(2.51)withrespecttotime,wesolvefortheradiusasafunctionoftimeR(t)=qR2022t:(2.52)Thelocationwhereisplacedinequation(2.46)fromotherreferences[14,44,54].Therefore,wepointoutthatourtimescaleshavebeenappropriatelyrescaledforcomparison.Themovinginterfaceproblemwassimulatedusing=0:05,t=6:41042=0:0256,andx=y=28ˇ0:0039,whicharebasedontheparametersusedin[44].ThenumericalsolutionisdisplayedinFigure2.4,andweobservethattheinterfacialcircleshrinks,asisexpected.InFigure2.5,weplotcomparetheevolutionoftheradiiobtainedbyoursecondorderschemewiththeexactradius(2.52),fortwotvaluesoftheparameter.Theradiusismeasuredbytakingasliceofthesolutionalongy=0,andthensolvingforthespatialpointwhereu=0usinglinearinterpolationbetweenthetwoclosestpointsthat41(a)u(x;y;0)(b)u(x;y;T2)(c)u(x;y;T)Figure2.4:Timeevolutionofthe2DACnumericalsolutionwithinitialcondition(2.50),uptotimeT=0:02562=10:24.(parameters:t=0:0256;x=y=0:0039;=0:05)satisfyu(x;0;t)<0andu(x;0;t)>0.ntisperformedwithaspatialmeshx=y=28,andtimestepsoft=0:2560;0:1280,and0:0640.Becausetheradiusisderivedasanexactsolutioninthelimit(i.e.,!0)[3],weobservethatthesmallervalueofisindeedmoreaccurate.WenextperformantstudyfortheAllen-Cahnequations,butthistimeintwospatialdimensions.Todoso,wemustincorporatethemultivariatesuccessiveconvolutionalgorithmsin(2.40)and(2.25)intothesecond(2.47)andthird(2.48)orderschemes.Giventhatwedonothaveanexactsolution,wecomputesuccessiveerrorsinanL1-norm.Thatis,wecomputejjutut2jj1foreachtimestept.Resultsareasexpected,andarepresentedinTable2.5.Theparametersusedare=0:05,x=y=29,andthecomputationtimeisT=0:5.Again,thequadraturemethodisfourthorderaccurateinspace,sothatthedominantsourceoferroristemporal.2.5.2FitzHugh-NagumoSystemFinally,wesolveawellknownreactionsystemthatarisesinthemodelingofneurons,theFitzhugh-Nagumo(FHN)model[28,39].TheFHNsystemconsistsofan42(a)=0:05(b)=0:01Figure2.5:Radiiofthenumericalinterfacialcircleasafunctionoftime(0t0:02562)withvarioustimestepsizes.(a)=0:05(b)=0:01(commonparameter:x=28).Redline:theexactradiusRin(2.52).P=1P=2P=3tL1errororderL1errororderL1errororder0:00636:59411041:17401045:17441060:00313:30651040:99593:26371051:84687:83511072:72340:00161:65631040:99738:67261061:91201:08111072:85740:00088:28941050:99872:23891061:95371:29611083:0602Table2.5:tstudiesfor2DAllen-Cahnequataionwiththeinitial(2.50)andwithhomogeneousNeumannboundaryconditions.activatoruandaninhibitorv,whicharecoupledvianonlinearreactionequationsut=Dur2u+1h(u;v);vt=Dvr2v+g(u;v);(2.53)whereDu,Dvarethecotsforuandv,respectively,and0<˝1isarealparameter.WeusetheclassicalcubicFHNlocaldynamics[39],thatareas43h(u;v)=Cu(1u)(ua)v;g(u;v)=udv;(2.54)whereC;aanddaredimensionlessparameters.Theparametersweusearethesameasin[19,50]:Du=1,Dv=0,a=0:1,C=1,d=0:5,and=0:005.ThecotfortheinhibitorisDv=0,identicaltotheworkfoundin[41,59].Thesecondorderschemefrom(2.47)isappliedtoeachvariableuandvseparately.Thisthenumericalschemeasun+1=etr2un+t2hn+t2hn+1;vn+1=vn+t2gn+t2gn+1;(2.55)wherehn=h(un;vn)andgn=g(un;vn).Weagainuseastabilizedpointiterationtoaddressthenonlinearreactionterms.Because(u;v)=(0;0)istheonlystableexcitablepointofequation(2.53)simplyconstructtheJacobianofF=(h;g)of(2.54)aboutthispoint:JF(u;v)2664uuvv3775@(h;g)@(u;v)j(0;0)2664uv3775=2664Ca11d37752664uv3775=2664Cauvudv3775:(2.56)44Theresultingsecondorderschemeis26641+Cat2t2t21+dt237752664un+1;k+1vn+1;k+13775=2664etr2un+t2hnvn+t2gn3775(2.57)+t226641hn+1;k+Caun+1;k+vn+1;kgn+1;kun+1;k+dvn+1;k3775;wherekistheiterationnumber,andnisthetimelevel.2.5.2.1Numericalresult:2DFitzHugh-NagumotestWesolvetheFitzHugh-Nagumosystemoverthedomain=[20;20][20;20]withPeriodicboundaryconditions.Usinginitialconditionsu(x;y;0)=8>><>>:0iffx<0g[fy>5g1(1+e4(jx5))2+1(1+e4(jx1))2otherwise(2.58)andv(x;y;0)=8><>:0:15iffx<1g[fy<10g0otherwise(2.59)wepresentthenumericalevolutionoftheactivatoruinFigure2.6.Weobservesimilarspiralwavesthatforminotherrecentworkfromtheliterature[19].45(a)T=1(b)T=2(c)T=4Figure2.6:Contourplot:Temporalevolutionoftheconcentrationofactivatoruusingthesecond-orderschemein(2.57)withparameters:t=0:01;x=y=0:2;M=4(fourthorderinspace).46Chapter3MOLTforCahn-HilliardEquation3.1ModelsInthisChapter,wepresenttheproposedmethodoflinestranspose(MOLT)schemeforCahn-Hilliard(CH)equationut=h2uf(u)i;x2ˆRd;(3.1)withtheappropriateinitialandboundaryconditions.Thephasefunctionu2H1describesthevolumefractionofonecomponentofabinarymixturewhereH1isthestandardsobolevspace.Thereactionfunctionf(u)isthederivativeofclassicalGinzburg-Landaudouble-wellpotential,F(u)=14(u21)2,whoselocalminimaisatu=1.Thesmallparameter0<˝1isthewidthoftheinterfacialtransitionlayer.Ina1Dsetting[a;b]ˆR),theLaplacianintheaboveequtionsisreplacedby@xx.In1893,vanderWaalsdevelopedalocalfreeenergydensityfunctionofbinarymix-tures,whichbecameakeycomponentinaldmodel.Overdecades,in1950s,CahnandHilliardrediscoveredthisfunctionsroleandextendedthisconcepttomodelspinodaldecompositionbasedontheiranalysis.Sincethen,thiscontinuumequation,Cahn-Hilliardnonlinearequation,hasbeenextensivelyusedtopredicttheevolutionofarbitrarymorphologiesandcomplexmicrostructureswithoutexplicitlytrackingthepositionsofinter-47faces[15].TheCHequation(3.1)describestheH1gradientwoftheCHfreeenergy[8]E(u)=Z22jruj2+F(u)dx:(3.2)Subjecttoboundarycondition,theCHequationdissipatestheCHfreeenergyddtE(u)=˝ut;Eu˛L20:Weemphasizethatthephysicalpropertyofenergystabilityisvital,andanyconsistentnumericalschemefortheCHmodelmustalsopossessthisproperty.Manysuchschemeshavebeenproposed,suchastheoperatorsplittingapproachintroducedbyEyre[26],thesemi-implicitspectraldeferredcorrectionmethodbyShen,[54,46],andevenfullyimplicitschemesutilizingnonlinearsolvers,suchastheconjugategradientmethod[16].WeshallusealinearlyimplicitpointmethodtosolvetheCHequation,whichisstabilizedbyshiftingthephaseaboutthebackgroundstateu=1.Inaddition,wealsoconsidervectorversionofCH(VCH)equation[16].Foru=(u1;u2),ut=h2uruW(u)i(3.3)wherethereactiontermisthederivativeofthepotentialfunction,W(u)=3i=1juzij2;zi:cuberootsofunityinthe(u1;u2)plane:(3.4)ThepotentialWisnon-negativeanditsminimumvaluesareattainedatthreevectorsfzig,48soastomodelathree-phasephysicalsystem.Thismodelcanbeseenasthehigherordervolumepreservingversionofthevector-valuedGinzber-Landauequation[5],whichsuggestsamodelforthreephaseboundarymotion,suchasthegrain-boundarymotioninalloys.Theenergyfunctionalisavectorversionof(3.2)suchthatEVCH(u)=Z22jruj2+W(u)dx:(3.5)3.1.1transformedmodelBeforeemployingatimediscretizationtoCHequation(3.1),weintroduceanewtrans-formedvariablev=u+1.Sincepurestatesu=1dominatesthesolutionduringripening,weseethatequivalentlyv=0;2.Thistransformationispartiallymotivatedby[18],inwhichthesameapproachiscombinedwithoperatorsplittingtoproduceacontractiveoperator.Similarly,in[16]thelinearizedvariablevˇu+1isusedasapreconditionerforthefullproblem.Ourworkisdistinct,howeverinthatwearenotsolvingthelinearizedequation,butinsteadthefullynonlinearequation,vt=2vxxxx+f(v1)xx;f(v1)=v33v2+2v;(3.6)whichfollowsfromthetransformationv=u+1insertedintotheCHequation(3.1),with=(a;b).Similarly,theenergyfunctional(3.2)becomesE(v)=Zba22jvxj2+F(v1)dx;F(v1)=14v4v3+v2(3.7)49InSection3.2,wewillformulateourMOLTschemeusingaboveequation(3.6),andwillshowthatthistransformationmakesourschemegradientstableinSection3.2.1.InthepreviousChapter2,wedevelopedMOLTschemeforHeatequationandAllen-Cahnequation.Wewillrevisitthispreviousworkwhennecessary.3.2Firstorderschemefor1DCHequationWeutilizetheMOLTbydiscretizing(3.6)intimeastheBackwardEuler(BE)scheme,vn+1vnt=2@xxxxvn+1+2@xxvn+1+@xx~fn+1;~f(v)=v33v2;(3.8)wherevn=u(x;tn)+1andt=tn+1tn.Thisscheme(3.8)isintimebutstillcontinuousinspace.Wenotethattheequationcontainsbothlinearandnonlinearimplicitterms,andsoantiterativeschemeisrequiredtoconstructafullyimplicitsolution.Wewillsuggesttwononlineariterationschemestosolvethesolutionvn+1inSection3.2.2.Wewouldliketogiveasimpleproofthatthesemi-discretesolutionof(3.8)isunconditionallygradientstablewithoutconsideringthespatialdiscretizationinthefollowingSection.3.2.1EnergystabilityInthisSection,wewillgiveasimpleproofthattheorderscheme(3.8)isgradientstableintheH1norminthesemi-analyticsettingunderthereasonableassumptions.Theproofisverysimilartothestabilizedsemi-implicitschemein[54]andoperatorsplittingschemein[26],exceptourschemeisfullyimplicit.Thesalientpointisthatonecanprovesuchstability50whentheoriginalsystemismothroughoursimplechangeofvariableinSection3.1.1.Wemakeusefulobservationswhichwillbeusedtoprovestability.Theoperator1:v2L2jRvdx=0!H2is[18]1v=h()<v;q>=<h;q>;8q2L2:(3.9)Foranyp;q2L2onecaneasilyprovefollowingidentity,<p;pq>=12kpk20kqk20+kpqk2012kpk20kqk20;(3.10)where<;>isstandardL2innerproductandkk0istheL2norm.Wewillconsiderthescheme(3.8)ingeneralspatialdimension2Rd;d2N).Lemma1.Underthefollowingassumption,0v2;v=u+1thefully-implicitscheme(3.8)thediscreteenergylawforanytimestept:E(vn+1)E(vn);8n0:wherevnisapproximationofv(x;tn)u(x;tn)+1ofCahn-Hilliardequation(3.1).Proof.Thesolutionvn+1thefollowingweakformulationinH21t<vn+1vn;˚>=2<vn+1;˚>+2<vn+1;˚>+<~fn+1;˚>;(3.11)51forall˚2H2Wechoose˚=1(vn+1vn),thenthisweakformbecomes1tkr1(vn+1vn)k20=2<rvn+1;r(vn+1vn)>+2<vn+1;vn+1vn>+<~fn+1;vn+1vn>;(3.12)bythenitionofoperator1in(3.9),andemployingintegrationbyparts.Wenowapplytheproperty(3.10),022(krvn+1k20krvnk20)+(kvn+1k20kvnk20+kvn+1vnk20)+<~fn+1;vn+1vn>:(3.13)WereplacethelasttermwithitsTaylorexpansion<~Fn+1~Fn;1>=<~fn+1;vn+1vn><~f0(˘n+1)2(vn+1vn);vn+1vn>;sothattheinequality(3.13)becomesE(vn+1)E(vn)=22(krvn+1k20krvnk20)+(kvn+1k20kvnk20)+<~Fn+1~Fn;1> k~f0k121!kvn+1vnk200;wherekk1isL1norm.Thelastinequalityholdsbecause~f0=3v26v0withourassumptionv.Therefore,withthisphysicalassumption,theimplicitscheme(3.8)isunconditionallygradientstable.Remark8.In[54,46],Shenet.al.introducedastabilizingtermS(un+1un)fortheirsemi-implicitschemefortheACandCHequations.IfSmaxu2Rj~f0(u)j2withthetruncated52potential,thentheirschemeguaranteesunconditionallystability.Inpractice,iftheyaretruncatingF(u)outsidetheinterval[1;1]suchthatmaxu2Rjf0(u)j=2.Essentially,ourofvisequivalenttothisassumption.Moreover,ourtransformedvariablehassamestabilizingct,withoutintroducinganadditionalterm.Weprovedtheweakformulationofourscheme.Inpractice,wesolveCHequa-tionusingastrongformulation.Eventhoughtheproofwegaveithereisonlyillustrative,allnumericalsimulationsweareusingallowstotakelargetimestepsizewhilemaintaininggradientstability.3.2.2TwononlineariterativeschemesInthisSection,wewillgivedirectiterationsolverstosolvethenonlinearproblem(3.8).Wesetupasimplepointiterationtosolvevn+1oftherscheme(3.8)It@xx2t@xxxx)vn+1;k+1=vn+t@xx~fn+1;k;(3.14)wherekindicatestheiterationindex.Bylaggingthethenonlinearterm~fn+1;k,theiterationupdateismadeexplicit.Now,weinvertthefourth-orderoperatoranalytically,yieldingtheupdateequationvn+1;k+1=It@xx+2t@xxxx1hvn+t@xx~fn+1;ki:(3.15)IfweemploytheGreen'sfunctionmethoddirectlytothisfourthorderoperator,thentheGreen'sfunctionwouldcontainbothdecayingandoscillatorycomputations,whichwouldnotproduceantmethod.However,wecaninsteadfactorthisintoaproductof53secondordertialoperators(modiHelmholtzequations)whichhaveanonoscilla-toryGreen'sfunction,respectively,andindoingsoleveragethefastcomputationmethodspreviouslyemployedinthepreviousChapter2.Aftersomealgebra,werewrite(3.15)astheproductoftwosecondordertialoperators(moedHelmholtz),vn+1;k+1= I@xx21!1 I@xx22!1hvn+t@xx~fn+1;ki;(3.16)where121=t+pt22tand122=tpt22t,sothat1;22Rwhent2.Fort2,theinverseoperatorsbothcorrespondtoGreen'sfunctionwhichdecayexponentially.However,fort<2,thecorrespondingGreen'sfunctionexhibitsoscillatorybehavior.Toovercomethisy,wesuggestanotherointiterationtosolvevn+1:vn+1;k+1=Ip2t@xx2"vn+t@xx fn+1;k2r2tvn+1;k!#;(3.17)whichisbasedoncompletingthesquare,andhasrequiredtheadditionoftheterm2p2tvn+1;k+1vn+1;k.ThepurposeforaddingtheextratermistoensurethattherequiredGreen'sfunctionmaintainsdecayingproperty.Weshouldpointoutthateventhoughtheiteration(3.16)hasatimestepsizere-strictiontohaveareal-valuedGreen'sfunction,wegaveasimpleproofthatthisschemeisgradientstableinsemi-analyticsettinginSection3.2.1.Moreover,wewillpresentnu-mericalevidencethattheschemeconvergesevenwithverylarget,whichisofparamountimportancefornumericalsimulationofCHtypeequations.54Figure3.1:Bifurcationoftheparameter12inthemoHelmholtzoperatorL.Whent2,alternatively,wecantransitionto(3.17).Sincebothiterationschemesconvergetoonesolutionwiththeerrorofsomeorderoft,andareexactlyequivalentatthetransitionpointt=2asshowninFigure3.1.Therefore,wesummarizeasfollowing:Factorizationmethot2):Invert I@xx21! I@xx22!Completedsquaremethodt<2):InvertI@xx22wheretheconstantsandi(i=1;2)aredeterminedbyeachschemeinFigure3.1.Weshouldpointoutthatthemainreasonwechoosetochangetheoperatorisbecausethesolutionforvisbetween0and2andanon-oscillatorykernelhelpsensurethattheintegralsolutionwecomputeisalsonon-oscillatory.Remark9.Theseconditerativescheme(3.17)wasdbyintroducingasecond-orderterm,S@xx(vn+1;k+1vn+1;k)whereS=2p2t.Sinceourvariablevwasalreadytrans-formedbythesteadystae,u=1,thisadditionaltermcouldlowerthestabilizingct.Indeed,wehavecheckedthatthisseconditerationguaranteesadecreaseinenergywhen55tpwherepˇ1:5,sothatwecansafelytransfertheiterationwithourcriteriat2.BothiterationsaremadebytwoinversionsofamoHelmholtzoperatorin(2.3)where=1and=2,henceL=I@xx2;withthentvaluesof.FormalinversionofthisoperatorinSection2.2.1yields:L1[v](x)=I[v](x)|{z}ParticularSolution+Bae(xa)+Bbe(bx)|{z}HomogeneousSolution;axb;wheretheparticularsolutionisaconvolutionwiththeGreen'sfunction,I[v](x)=2Zbaejxx0jv(x0)dx0;andthecotsBaandBbaredeterminedbyapplyingboundaryconditionsatx=a;binSection3.2.3.Remark10.Weadditionallycancalculatethesecondderivativeofuasfollows@xx=2(IL)=)L1[@xxv]=2L1I[v];whereweusetheofentialoperator(2.3).3.2.3FullydiscretesolutionOursemi-discreteschemeismadebydoubleinversionofmoHelmholtzoperators,inSection3.2.2.ItremainstodiscretizetheconvolutionintegralI[v](x),andobtainafullydiscretealgorithm.InthepreviousChapter2,wehaveaccomplishedthisinversion56withthefastconvolutionalgorithmforourinversionoperator(2.9);forconvenience,wewillpointouttheessenceofthisfastalgorithminthisSection.Forexample,wereviewthespatialquadratureonauniformgrid(hj=x),asthesecond-orderaccurate(M=2)scheme.Firstofall,thedomain[a;b]ispartitionedintoNsubdomains[xj1;xj]a=x0<x1<<xN=b;hj=xjxj1;wethenshouldevaluatetheconvolutionoperatorateachgridpointthroughIjI[v](xj)=ILj+IRjinSection2.2.2.1.However,becauseILandIRaresatisfyingtherecursivere-lationin(2.12)and(2.13)respectively,thisevaluationcanbethroughfollowingapproximationILj=eILj1+JLj;JLj=2Z10ezv(xjxz)dz;IRj=eIRj+1+JRj;JRj=2Z10ezv(xj+xz)dz;whereJLjandJRjare"local"integralsinSection2.2.2.2and=x.ThestandardLagrangeinterpolationisusedtoapproximatetheselocalintegrals,JLjˇZ10ezpLjdz;pLj=z2+z2vj1+1z2vj+z2z2vj+1;JRjˇZ10ezpRjdz;pRj=z2z2vj1+1z2vj+z2+z2vj+1:57Integratingagainstanexponential,wethatJLjˇw1vj1+w1vj+w2vj+1;JRjˇw1vj1+w0vj+w1vj+1;wherethequadratureweightsaregivenbyw1=12˚j;2+˚j;1=12d+13d2+1d2;w0=˚j;0˚j;2=12+d1d2;w1=12˚j;2˚j;1=121d21+d2:whered=eand˚j;kisin(2.17).Forhigher-orderspatialschemes(M>2),thequadratureweightsarepre-computedasthesimilarway,whichpresentedinSection2.2.2.2,sothatthefastconvolutionalgorithmisachievedasO(MN)pertimestepwhereorderMinspace.Ineverysimulationinthiswork,wechooseM=4orM=6.Next,thehomogenoussolution(2.11)isusedtoenforcevariousboundaryconditionsinSection2.2.3.3.Forexample,periodicboundaryconditionsleadtovn(a)=vn(b);vnx(a)=vnx(b);8n2N;andevaluationof(2.11)atx=a;b,producesa22systemfortheunknownBaandBb.SolvingthislinearsystemyieldsBa=IN1;Bb=I0158where=e(ba)andI0=I[v](a)andIN=I[v](b).3.3Higherorderschemefor1DCHequationInthisSection,weshowhowtomakesecondandthirdordertimeaccuratemethodsbycombiningtheideasofMOLTformulationwithBackwardFormulas(BDF),SinglyDiagonalImplicitRungeKutta(SDIRK),andSpectralCorrection(SDC)methods.Wewillpresentseveraltstudiestocomparethosemethods.3.3.1MOLTwithBackwardFormula(BDF)TheBDFtimesteppingmethodsareoneofthemostcommonlyusedimplicitLinearMulti-stepMethods[34].Toachievethedesignedorderofaccuracy,wediscretizeatimederivativeusingBDFformulasasfollows:BDF2:I+232t@xxxxvn+1=43vn13vn1+23t@xxfn+1;(3.18)BDF3:I+6112t@xxxxvn+1=1811vn911vn1+211vn2+611t@xxfn+1;(3.19)wherewenowrequiretwoorthreeprevioustimesteps,respectively.SimilarwiththeorderschemeinSection3.2,wetwoointiterations.First,energy-stablefac-torizationiterationwhichhasatimesteplowerboundforreal-valuedGreen'sfunctions.Second,completedsquareversionbyaddingS(vn+1;k+1vn+1;k)term,whichhasatime59stepupperboundforenergystability.Withthat,BDF2iterationsfor(3.18)t322:L=I43t@xx+232t@xxxx= I@xx21! I@xx22!;t<322:L=I2r232t@xx+232t@xxxx=I@xx22;where121=23t+q(23t)2232t,122=23tq(23t)2232tand12=q232t.ForBDF3in(3.19),t1162:L=I126t@xx+6112t@xxxx= I@xx21! I@xx22!;t<1162:L=I2r6112t@xx+6112t@xxxx=I@xx22;where121=611t+q(611t)26112t,122=611tq(611t)26112tand12=q6112t.ThemainadvantageofBDFschemesisthattheextensionofmethodstohigher-orderonesisverystraightforward.However,theyalsorequireinitializationofthefewtimesteps.InnextSection,wewilllookattheRugeKutta(RK)method,whichdonothavethisrequirement.3.3.2MOLTwithSinglyDiagonalImplicitRungeKutta(SDIRK)TheSDIRKmethodisanimplicitRungeKutta(RK)method,whichhasthesamediago-nalelementintheButcherTablein[2].ForthePthordermethod,RKmethodsneedPintermediatestepsinonetimestepping.Itissomewhatcomputationallyexpensive,butitisdesirablefortheinitialstepofmultistepmethods.Here,wederivethesecond-orderSDIRK60(SDIRK2)methodforCHequation,butemployingSDIRK3issimilar.I+2t@xxxxK1=2@xxxxvn+@xxf(vn+tK1);=1p22;I+2t@xxxxK2=2@xxxx(vn+(1tK1)+@xxf(vn+(1tK1+tK2)vn+1=vn+tf(1)K1+K2g;whereistherooltofthepolynomial122+2,whichisderivenfromtheorderconditions[2].Next,wesolvetwointermediatesolutionK1andK2usingsimilarnonlineariterativeschemesinSection3.2.2.First,ift2,Kn+1;k+11=L11L12h2@xxxxvn+@xx~fvn+tKn+1;k1+2vn;Li=I@xx2iKn+1;k+12=L11L122@xxxxvn+(1tKn+1;k1+@xx~fvn+(1tKn+1;k1+tKn+1;k2+2(vn+(1tKn+1;k1);where121=t+p(t)22tand122=tp(t)22t.Ift<2,Kn+1;k+11=L22@xxxxvn+@xxfvn+tKn+1;k12q2tKn+1;k1;Kn+1;k+12=L22@xxxxvn+(1tKn+1;k1+@xxfvn+(1tKn+1;k1+tKn+1;k22q2tKn+1;k2;whereLIp2t@xxisused.61Remark11.SimilartoRemark10,thefourthderivativecanbecalculatedasfollows,@xxxx=2122(IL1)(IL2)=)L11L12[@xxxxv]=2122L11IL12I[v]:InSection3.3.4,wewillpresenttstudiesforboththeBDFandSDIRKmeth-ods.Wenotethatforsmallenoughtimestept,bothmethodsconvergeasexpected.However,forthecorrespondingthird-ordermethods(BDF3andSDIRK3),astincreasestheorderofconvergencebeginstoplateau.Toresolvethisissue,wewillsuggestanotherhigherordermethod.3.3.3MOLTwithSpectralDeferredCorrection(SDC)TheSDCmethodsareaclassoftimeintegrators[23].DeferredCorrection(DC)methodscomputeapredictiontothesolution("level0")usingloworderschemes(e.g.BackwardEuler),andthencomputeoneormorecorrectionsatthesubsequentlevels,togethigherorderofaccuracy.WepresentthetraditionalSDCprocedurein[17].1.Predictionstep(level[0]):Subdividethetimeinterval[0;T]withuniformt:0t0<t1<<tNtT:Computefv[0]m(x)g0mNt,viaBackwardEulerapproximationin(3.8)forCHequa-tion,whichisapproximationtotheexactsolutionfy(x;tm)ym(x)g0mNt.622.Correctionstep(level[1]):Assumethatv(0)(x;t)beapolynomialinterpolationtotheexactsolutiony(x;t)satisfyingv(0)(x;tm)v[0]m(x);m=n+1;n;;n+1P;wherePwillbesplater.Theerrorequationisde(x;t)=y(x;t)v(0)(x;t);andtheresidual(or"defect")is(x;t)=v(0)tFCH(v(0));FCH(v)=2@xxxxv+@xxf(v):Wetakethederivativeoftheerrorequationwithrespecttot,andrewriteitusingtheresidualet(x;t)=yt(x;t)v(0)t(x;t)FCH(y(x;t))FCH(v(0)(x;t))(x;t);()et(x;t)+(x;t)=FCH((e+v(0))(x;t))FCH(v(0)(x;t));()@@tv(1)(x;t)Zt0FCH(v(0)(x;˝))d˝=FCH(v(1)(x;t))FCH(v(0)(x;t));whereweassumethattheinitialconditione(x;0)=0,andv(1)=v(0)+e.Hence,ourupdating"level[1]"solutionfv[1]mgisfoundbyapproximatingtheabovetialequationasthesamemethodas"level0".Inourcase,weapplyBackwardEuler,v[1]n+1v[1]n=tFCH(v[1]n+1)FCH(v[0]n+1)+Ztn+1tnFCH(v(0)(x;˝))d˝:(3.20)633.Correctionstep(level[j]):Theprocessistheniteratedbygeneralizing(3.20),v[j]n+1tFCH(v[j]n+1)=v[j]ntFCH(v[j1]n+1)+Ztn+1tnFCH(v(j1)(x;˝))d˝;(3.21)wherethetermsincludingupdatingsolutionv[j]havebeencollectedonthelefthandside,andtheoldsolutionv[j1]areontherighthandside.TocompleteSDCmethod,wemustconsideranapproximationoftheintegralin(3.21):Ztn+1tnFCH(v(j1)(x;˝))d˝=8>>>>><>>>>>:tPXi=0~qiF[j1]n+1i;ifnP1;tPXi=0~qiF[j1]iifn<P1;(3.22)whereF[j1]n=FCH(v[j1]n),and~qiarequadratureweights(cf.[17]).Forexample,ifP=2andn>P1,thefunctionvaluesF[1]n+1,F[1]n,andF[1]n1areusedforthequadratureand~qi=Ztn+1tnP=2Yk=0;k6=i(ttn+1k)(tn+1itn+1k)dt;i=0;1;2:Notethatthenumberoftermsinthesum(3.22)isP+1,wherePistheorderofpolynomialinterpolationv(j1).Theintegralmustbeapproximatedwithincreasingaccuracyaslevelincreases,sothatPjat"level[j]".InSection3.3.4,wewillshowthatthisorderPtheasymptoticregionofstabilityofSDC3method.64WenowcombineMOLTschemewiththesehigher-orderSDCmethods.Forinstance,thesecond-orderSDCscheme(SDC2)onlyrequiresonemorecorrection(level[1]).Ift2,v[1]n+1;k+1=L11L1224v[1]n+t@xx0@~f[1]n+1;kf[0]n+1f[0]n21A+2t2@xxxxv[0]n+1v[0]n35(3.23)wherequadratureweightsof(3.22)are~q1=~q2=12(P=1:trapezoidalrule)andLiaresamewiththeschemeinSection3.2.Similarly,ift<2;v[1]n+1;k+1=L224v[1]n+t@xx0@f[1]n+1;k2r2tv[1]n+1;kf[0]n+1f[0]n21A35+L22t2@xxxxv[0]n+1v[0]n(3.24)Wewillnowpresenttstudiesforsecond-andthird-orderBDF,SDIRKandSDCmethodsinthefollowingSection.3.3.4Numericaltest:tstudiesof1DCHsolutionsInthisSection,wewillchecktheorderofaccuracyofpresentedtimesteppingmethodsandcomparetheseschemes.Startingfromaninitialdatainreference[16],u0(x)=cos(2x)+1100ecosx+110;x2[0;2ˇ](3.25)andwiththeperiodicboundaryconditions,weintegrateuptoatimeTusingeachsecond-ordermethod.Ifwedenoteutasthecomputedsolutionwithtime-stepsizet,thenrecomputesolutionut=2withhalvedtimestepuptoT.Giventhatwedo65BDF2SDIRK2SDC2tL1errororderL1errororderL1errororder0:05002:94541054:77871066:37631050:02507:27711062:01701:20501061:98761:54121052:04870:01251:97551061:88123:00811072:00213:76181062:03450:00635:38691071:87477:50761082:00249:27861072:01940:00311:43521071:90821:87321082:00292:30321072:0103Table3.1:tstudiesofsecond-ordermethodsfor1DCHequationwithperiodicBC.nothaveanexactsolutionof(3.1),weestimatedtheerrorfortasthesuccessiveerrorinmaximumnorm,thatiskutut=2k1,andagainhalvedthetimesteptocomputethecorrespondingerrorrepeatedly.Thosesuccessiveerrorsforthreesecond-ordermethodsarepresentedinTable3.1.Weusedtheparametersincommon:=0:18;x=2ˇ512ˇ0:0123;T=1;Ntol=1012;(3.26)whereNtolisatoleranceforpointiterationssuchthatkvn+1;k+1vn+1;kk1<Ntol.Weuse6th-orderspatialquadrature(M=6),sothatthedominanterrorwastemporal.Asexplained,eachmethodhasacriterionfortimestepsize,whichleadstoswitchiterativeschemes:BDF2:t=322;SDIRK2:t=21p22;SDC2:t=2:Withthiscriteria,weswitchedthecorrespondingointiterations,andtheTable3.1showsthateachofthemethodsexhibitssecondorderconvergence,respectively.Thisresultguaranteesthesafeadaptiveswitchingcontrol.MoreoverinFigure3.2a,weseethatthetotalenergyofeachtimesteppingmethoddecaysduringthetimeevolution.66(a)Energydecay(2ndordermethods)(b)Iterationcount(2ndordermethods)(c)Energydecay(3rdordermethods)(d)Iterationcount(3rdordermethods)Figure3.2:Energydescentandnonlineariterationcounthistoryofallsecond-andthird-ordermethodswithtimestept=0:025andparameters(3.26).67InFigure3.2b,weplottheiterationcountateachtimelevelofallsecond-orderschemes(withthesametoleranceNtol=1012).SinceSDIRK2iscomprisedottwointermediatesteps,thisrequirestwoseparatenonlineariterations.Hence,SDIRK2methodsismoreex-pansivethantheothertwomethods.Infact,usingourMATLABcodes,thecomputationaltimeofBDF2,SDC2,SDIRK2is0:9725(s),1:9003(s),and2:6468(s),respectively,togetnaltimeT=1.Therefore,wecanconcludethatBDF2ismostcomputationallytamongtheseschemes.Next,wesimilarlyimplementedthecorrespondingthird-orderschemesforCHequation,withthesameparametersin(3.26),andwiththefollowingconsiderations:[BDF3]Needtwoinitialstep:SDIRK2used.(Switchcriterion:tswitch=1162)[SDIRK3]NeedtocomputeK1;K2;K3ateachupdate.tswitch=2;ˇ0:4359)[SDC3]Computefromlevel[0]tolevel[2]withthequadratureorderPin(3.22):IfwechooseP=2atlevel[2],thenweightsare~q0=512;~q1=812;~q2=112.tswitch=2)TheenergyanditerationhisotryofthirdorderschemesareshowninFigure3.2candFigure3.2d.Againweseethatthediscreteenergyofeachsolutiondecays.TheverybeginningiterationnumberofBDF3islargerthanBDF2'ssinceweinitializedusingSDIRK2.ButBDF3isstillthemosttmethodthanotherthirdorderschemes.UsingMATLABbuilt-infunction,thecomputationaltimeofBDF3,SDC3,SDIRK3is1:0355(s),2:8539(s),and3:6613(s),respectively,tocomputeuptothenaltimeT=1.Thesuccessiveerrorsofthethird-ordermethodsareshowninTable3.2.Eachmethodachievesthirdorderconvergenceforsmallt;buttheorderofconvergencebeginstoplateauforlargertimesteps.However,wepointoutthattheplateauislesssevere68BDF3SDIRK3SDC3tL1errororderL1errororderL1errororder0:05001:60011051:72591061:26691050:02505:28321061:59873:61181072:25662:10021062:59270:01251:17561062:16806:28231082:52343:11241072:75450:00631:93141072:60579:61631092:70774:26151082:86860:00312:65061082:86521:30111092:88585:41931092:97520:00163:25111093:02731:744710102:89866:691310103:0177Table3.2:tstudiesofthird-ordermethodsforthe1DCHequationwithperiodicBC.SDC3(P=3)tL1errororder0:05009:94621060:02501:30611062:92890:01251:41421073:20720:00631:61561083:12980:00311:78521093:17790:00162:275310102:9720Table3.3:tstudyofSDC3(P=3)forthe1DCHequationwithperiodicBC.forSDC3thanBDF3orSDIRK3;theconvergenceorderCR:=log2kutut=2k1kut=2ut=4k1ofSDC3startsaround2:59,whichisclosertotheexpectedorder,comparedtoBDF3(CR=1:6)andSDIRK3(CR=2:26)whent=0:05.Furthermore,ifweusehigherquadratureoforderPin(3.22)forSDC3,thentheantic-ipatedrateofconvergenceisobservedevenforlargetimestepst.Forinstance,ifweuseP=3,(thequadratureweightsare~q0=924;~q1=1924;~q2=524and~q3=124),thentherateisroughlythirdorderevenwiththelargertimesteps,showninTable3.3.WesummarizethisSectionthatwecanachieveexpectedorderofaccuracythrough(short-time)mentstudiesofsecondandthirdorderschemes.Fromthepointofviewofsolver,higher-orderBDFschemes(BDF2andBDF3)areideal.High-order69quadratureinSDC3showsinterestingresult,whichperformswellonlargertimesteps.3.4MOLTforMulti-DCHequationWenowextendthe1Dsolvertomultiplespatialdimensionviadimensionalsplitting[13,10,12]whichalsodiscussedinSection2.4.EachtimesteppingmethodpresentedthusfarrequirestheinversionofamoHelmholtzoperator,withindependently.Forsimplicity,in2D,L=I2=I@xx2I@yy2@xx@yy4LxLy@xx@yy4;(3.27)whereLxandLyareourO(N)1DsolversandthesubscriptsdenotethespatialcomponentofunivariatemoHelmholtzoperators.Later,weshouldformallyinvertbothoperators.NowL1xisdonefory,andviceversaforL1yinaline-by-linefashion,similartoalternatingdirectionimplicit(ADI)typemethod[22,27].However,thisfactorization(3.27)yieldsamixedderivativeterm,whichisO(14).Ourkeyobservationisthatifweincludethisterminourpointiteration,thenwecansimultaneouslysolvethenonlinearproblemandremovethesplittingerror.Moredetaildanalysisofthe2DCHequation(3.1)and2DCHvectorequation(3.4)willbepresentinnextseveralSections.3.4.1Firstorderschemefor2DCHequationInthisSection,wewillpresentthesemi-discreteschemefor2DCHequation,usingtheBackwardEulermethod.Theanalogoushigher-orderschemesfollowaccordingly.Startingfromthe1Dscheme(3.14),wereplace@xxwiththeLaplacianoperator=@xx+@yy.If70t>2,wefactorizetheoperator,andhave I21! I22![vn+1;k+1]=vn+t~fn+1;k;~f=f(v)2vv33v2;(3.28)where12i=tpt22t,i=1;2respectively,whicharesameparameterswith1Dscheme.Pluggingtheidentityin(3.27),bylaggingthemixedderivativetermalongwiththenonlinearterm,L1;xL1;yL2;xL2;y[vn+1;k+1]=vn+t~fn+1;k+ L1;xL1;y@xx@yy42+L2;xL2;y@xx@yy41 @xx@yy41! @xx@yy42!!vn+1;k;(3.29)whereLi;x=I@xx2i;Li;y=I@yy2i;(i=1;2).NotethatLaplacianoperatorandmixedderivativecanbereplacedasfollowing:=2ILxLy+@xx@yy4;@xx@yy4=(LxI)(LyI):(3.30)NowweformallyinvertbothoperatorsL1;xL1;yL2;xL2;ytotherighthandsideof(3.29),vn+1;k+1=L2;xL2;yL1;xL1;y1[vn]+L2;xL2;y1C1[21t~fn+1;k]+M1;2hvn+1;ki(3.31)71whereCi=Li;xLi;y1I+Di;xDi;y;D=IL1;(=fx;yg);(3.32)Mi;j=Di;xDi;y+Dj;xDj;yDi;xDi;yDj;xDj;y:(3.33)Asshown,everymixed-derivativesplittingerrortermcanbecontrolledbyapplyingD,(=fx;yg)operators,whichcanalsobeconstructedinaline-by-linefashion.WeemphasizethatthisallowsustoremovesplittingerrorO(14)=Ot).Similaroperators'extensionL1andDholdsfortheothercase,t2: I20!2[vn+1;k+1]=vn+ tfn+1;k220vn+1;k!;(3.34)where120=p2t.Thiscompletedsquareformisreplacedbytheidentityin(3.27),L0;xL0;y2[vn+1;k+1]=vn+ tfn+1;k220vn+1;k!+0@2L0;xL0;y@xx@yy40 @xx@yy40!21Avn+1;k:(3.35)Finally,weinvertthemoHelmholtzoperatorsandusetheaboveoperators(3.32)and(3.33),thenvn+1;k+1=L0;xL0;y2[vn]+L0;xL0;y1C0h20tfn+1;k2vn+1;ki+M0;0hvn+1;ki:(3.36)Wenowhaveswitching-schemewiththetime-stepcriteriont=2,whcihissameformulationwiththe1DinSection3.2.2.Inaddition,wecanextendthefully-72discretesolutionbyapplyingthesameprocedureinSection3.2.3vialine-by-linefashion.Thisstraightforwardextensionisoneofmainadvantagesofdimensional-splitalgorithm.WealsopointoutthatwehavealreadyprovedenergystabilityofaboverschemeinSection3.2.1.3.4.2Higherorderschemefor2DCHequationThehigher-ordertimesteppingmethods'extensionson2DareanalogoustotheorderschemeinpreviousSection3.4.1.Hence,wewillbstatetheformulasforsecond-order2DCHschemesinthisSection.BDF2 120=r232t;12i=23tr49t2232t(i=1;2)!t322:vn+1;k+1=L0;xL0;y243vn13vn1+L0;xL0;y1C02320tfn+1;k2vn+1;k+M0;0hvn+1;ki:t>322:vn+1;k+1=L2;xL2;yL1;xL1;y143vn13vn1+L2;xL2;y1C12321t~fn+1;k+M1;2hvn+1;ki:SDIRK2 120=p2t;12i=tp(t)22t;=1p22!73t2:Kn+1;k+11=1tC20[vn]+M0;0hKn+1;k1i+L0;xL0;y1C0h20f(vn+tKn+1;k1)2Kn+1;k1i:Kn+1;k+12=1tC20hvn+(1tKn+1;k1i+M0;0hKn+1;k2i+L0;xL0;y1C0h20f(vn+(1tKn+1;k1+tKn+1;k2)2Kn+1;k2i:vn+1;k+1=vn+(1tKn+1;k+11+tKn+1;k+12:t>2:Kn+1;k+11=1tC1C2[vn]+M1;2hKn+1;k1i+21L2;xL2;y1C1h~f(vn+tKn+1;k1)+2vni:Kn+1;k+12=1tC21hvn+(1tKn+1;k1i+M1;2hKn+1;k2i+21L2;xL2;y1C1h~f(vn+(1tKn+1;k1+tKn+1;k2)i+21L2;xL2;y1C1h2(vn+(1tKn+1;k1)i:vn+1;k+1=vn+(1tKn+1;k+11+tKn+1;k+12:SDC2:Thenumericalsolutionatpredictionstep(level[0])fv[0]ngn2NisobtainedbytheBEschemesinSection3.4.1.Thefollowingisthecorrectionstep(level[1])fv[1]ngdependingonthetimestepsize. 120=p2t;12i=tpt22t(i=1;2)!74t2:v[1]n+1;k+1=L0;xL0;y2hv[1]ni+L0;xL0;y1C0h20tf[1]n+1;k2v[1]n+1;ki+M0;0hv[1]n+1;ki12C20hv[0]n+1v[0]ni+20t2L0;xL0;y1C0hf[0]n+1f[0]ni:t>2:v[1]n+1;k+1=L2;xL2;yL1;xL1;y1hv[1]ni+L2;xL2;y1C1h21t~f[1]n+1;ki+M1;2hv[1]n+1;ki12C1C2hv[0]n+1v[0]ni+21t2L2;xL2;y1C1hf[0]n+1f[0]ni:Withtheseformulas,wewillimplementtstudiesfor2DCHequationinthenextSection.3.4.3Numericaltest:tstudiesof2DCHsolutionsInthisSection,wepresentresultsforthe2DCHequation.Weconsiderthestandardbenchmarkinitialstatesin[16,57]tothetemporalorderofaccuracyin2Dsetting.u0(x;y)=2e(sin(x)+sin(y)2)+2:2e(sin(x)sin(y)2)1;(x;y)2[0;2ˇ]2(3.37)withtheperiodicboundaryconditionsandwiththefollowingparameters,=0:18;x=2ˇ128ˇ0:0491;T=1(0tT);Ntol=106;(3.38)wherea4thorderspatialquadrature(M=4)isused.Wepresentatemporaltstudyofeachsecond-orderschemeispresentedinTable3.4.Thiststudyshowsthesecond-orderofconvergenceforallthreemethods.Furthermore,thetotalenergyofeachsolutionisdecreasedduringtimeevolution,asshowninFigure3.3a.Thedimensionalsplitalgorithmguaranteestheenergy-decentsolutioninmultiplespatialdimensions.BycomparingtheiterationcountateachtimelevelinFigure753.3b,theBDF2methodisthemostt,aswasthecasein1D.Wenotethatweaddressedthesplittingerrordirectlyinthedimensionalsplittingfor-mulationbyexplicitlyincorporatingthemixed-derivativeterm(viaDoperatorsin(3.32))intotheointmethodusedtosolvethenonlinearproblem.Thisnumericalstudiesguarenteethatwecanachievehigherordersofaccuracyinspaceandtimebyremovingsplittingerrors.Thisisimportantfeatureofthepresenteddimensionalsplitalgorithm.WewillshowthatoursolvercanbeeasilyappliedtothevectormodelinnextSection.3.5MOLTforvectorCH(VCH)equationWenowextendourdimensionalsplitMOLTalgorithmtovectorCH(VCH)model(3.3).Thismodelconsistsoftwocoupledvariablesu1andu2withlocaldynamicruW,comprisedofpartialderivativesof6thorderpolynomialsW(u1;u2)in(3.4).BeforeweemployMOLTformulationtothissystem,wetransformthevector(u1;u2)aboutitsbackgroundstate.Todothat,weneedtotheequilibriumpointsofthefunctionalW(u1;u2),(u1;u2)=f(cos(i);sin(i))j1=0;2=2ˇ3;=2ˇ3gf(1;0);(12;p32);(12;p32)gwhicharecuberootsofunityin(u1;u2)plane.AstraightfowardcalculationyieldstheJacobianofthepotentialruWatthosepointsJruW(u1;u2)26664@2W@u21@2W@u1u2@2W@u2u1@2W@u2237775(u1;u2)=26641800183775;(3.39)76BDF2SDIRK2SDC2tL1errororderL1errororderL1errororder0:10003:78911031:02501032:33341030:05008:26261042:19722:69501041:92725:05701042:20610:02502:03191042:02386:97961051:94911:10001042:20070:01254:69091052:11491:78091051:97062:83121051:9581Table3.4:tstudiesofsecond-ordermethodsforthe2DCHequationwithperiodicBC.(a)Energydecay(2ndordermethods)(b)Iterationcount(2ndordermethods)Figure3.3:Energydescentandnonlineariterationcounthistoryofallsecond-order2Dmethodswithtimestept=0:05andwithparameters(3.38).77thus,allpointsarestableequilibriumsolutions[57].Wesubractu=(u1;u2)oftheback-groundstatez3(12;p32),andintroducethenewtransformedvectorv=(v1;v2)=uz3intotheoriginalsystem(3.3),vt=22v+rvW(v+z3)22v+rv~W(v)+18v(3.40)where~Wv1(v):=Wv1(v+z3)18v1and~Wv2(v):=Wv2(v+z3)18v2.Wewilldevelopthetimemarchingschemeusingthistransformedsystem(3.40)inthefollowingSection.3.5.1TwononlineariterativeschemesWeapplytheBackwardEulerschemetothetransformedsystem(3.40),thenIt+2t2vn+1;k+1=vn+rv~Wn+1;k(3.41)wherekisaniterationindex.Weagainintroducetwofactorizationsoftheleft-handsideoperator:t281:It+2t2= I21! I22!;12i=tpt22t;(3.42)t<281:I2p2t+2t2=I22;12=p2t:(3.43)WecanhandleaboveoperatorswiththesamestrategyinSection3.4.1.ThehigherorderexpansionsareanalogousinSection3.4.2toeachcomponentv1andv2respectively.Asexpected,thestabilizededpointiteration(3.42)alsoallowustotakelargetimesteps78withoutforgoingtheenergystabilityproperty(3.5).3.5.2Numericaltest:tstudiesof2DVCHsolutionsInthisSection,weimplementoursecond-orderschemeto2Dsystem(3.3)forvectorCHmodel.Inordertolookspinodalevolutionofphasefunctionu=(u1;u2),webeginwiththesameinitialcondition(3.37)foru1andu2(x;y;0)=sin(y);(x;y)2[0;2ˇ]2(3.44)foru2inthereference[16].Withtheparameters=0:32;x=y=2ˇ128ˇ0:0491;Nmaxit=500;Ntol=106;(3.45)weimplementBDF2,SDC2,andSDIRK2methodsforVCHsystemduring0tT=0:1withtimestepsizet=0:000125.(a)Energydecay(2ndordermethods)(b)Iterationcount(2ndordermethods)Figure3.4:Energydescentandnonlineariterationcounthistoryofallsecond-orderVCH2Dmethodswithtimestept=0:000125andwithparameters(3.45).79AsshowninFigure3.4a,everyschemepreservesthediscreteformoftheenergylawduringthespinodalevolutionofvectorCHmodel.Figure3.4bisobtainedbyiterationcountateachtimesuchthatmaxkun+1;k+11un+1;k1k1;kun+1;k+12un+1;k2k1<Ntol.TheresultindicatesthatBDF2ismostt,asexpected,thusitrequiresleastnumberofnonlineariterationwhileallhavesameconvergencetoleranceNtol=106.Infact,computationaltimeusingbuilt-inMATLABfunction(tic-toc)ofBDF2simulationwas161.09seconds,whileSDC2is373.74secondandSDIRK2is803.15seconds.Inthissimulation,weusedsmalltimestepsizetoaccratelysimulatethespinodalphase.However,wealsopresentnumericallythatourschemesallowlargetimestepsanditerativeschemesconvergeinSection3.7.3.3.6TimeadaptivestrategyFormodels,adaptivetimesteppingisacrucialfeatureforantandaccuratenumericalsolution.Forinstance,thesolutionofCHequation(3.1)evolvesonvarioustimescales.Duringspinodalevolution,transitionlayersaredevelopedinO(1)time.Subsequently,theyslowlyevolveandmergeonalongertimescale,O(eC)for1Dmodel,whichiscalledripeningprocess.Simulatingthesephenomenawithatimesteppingnecessarilybecomet,andsowesuggestanadaptivetimesteppingstrategycombinedwithpresentedMOLTschemes.TheadaptivetimestepsizecontrolisbasedontheLocalTruncationError(LTE)eateachtimelevelt=tn.TheLTEcanbeapproximatedbyeˇ=jjuunjj1,whereuisanexplicitpredictorsolution,typicallyusingtheForwardEuler(FE)orAdamsBashforth(AB)schemesin[16].Ontheotherhand,theRichardsonextrapolation(known80Algorithm1Adaptivetimestep-sizecontrol1.Startingatt=tn,approximatethelocaltruncationerror.Forinstance,forRichard-sonextrapolation,solutionisestimatedtwice:onceasafullstepwitht(denotedbyun+1t,andastwohalfsteps(un+1t=2).ThencebetweenthetwonumericalapproximationsgiveanestimateforLTEofun+1(e.g.Richardsonextrapolation)eˇ:=12jjun+1tun+1t=2jj1:2.eatolerance˙tolfortheaboveLTE.Ifaccuracyfails(>˙tol),thenthetimestepfailsandrepeatwiththereducedtimestep.Ifaccuracysuccess,thenwealsotestthefollowingcriteriafortimestep-sizeselection,(I)˙tol:NitNmaxit<0:7:tn+1=tnminr˙tol;;=0:8;=1:3>10:7NitNmaxit<1:tn+1=tnminr˙tol;1;NitNmaxit1:Stepfails.Trywithshrinktn=tn1(II)>˙tol:Stepfails.Trywithshrinkstep-size˙tol>2:tn=tn1;˙tol2:tn=tnr˙tolwhereandaresafetyfactorsin[16].81asstep-doubling)orembeddedRunge-Kuttapairscanbeused[17].WeadoptRichardsonextrpolationaswelloursecond-orderscheme,BDF2andSDC2,toapproximatetheLTE.Wealsousethetimestep-sizeselectioncriteriapresentedin[16],whichissummarizedbelowAlgorithm1.Inpractice,theprocedureleadstosmalltimestepsduringspinodalevolution,orattheripeningevent,tomaintainaconsistentLTE.Ontheotherhand,duringslowcoarsening(metastablestates),smalltimestepsareunnecessaryandsotisincreasedwithintheupperboundforointiterationcount.3.7NumericaltestsInthisSection,wepresentadaptivetimesteppingresultsusingthepreviouslydevelopeddimensional-splitMOLTschemes.AsmentionedinSection1.2.1.3,wewillreproducenovelbenchmarkproblemsin[16]usingourschemes.Wewillcompareseveralxedtime-steppingschemeswiththeadaptivetimesteppingschemesvianumericalresults.3.7.11DCahn-HilliardModelInthisSection,wesolvethe1DCHequation(3.1)withainitialcondition(3.25)(=0:18).Thesecondperturbationtermofthisinitialstatecreatestwointervals,u=1andu=+1,areasymmetric,sothatenumberoftransitionlayersareformedduringspinodalevolution.Afteralongripeningprocess,suchlayersareeventuallyabsorbedintooneregion,attheso-calledripeningtime[16].Theaimofthissimulationistoaccuratelycapturealltimescalesusingbothandadaptivetimesteppingstrategies.(Thespatialstep-sizeiswithx=2ˇ128ˇ0:05.)82Intheexperiment,weimplementourvarioustimesteppingschemes,withsmalltimestept=0:01).TheripeningtimeTrisedasthatforwhichthemidpointvalueu(ˇ;t)changesfrompositivetonegative.ThepointiterationhasresidualtoleranceNtol=1011ateachstep,andtheripeningtimesarepresentedinTable3.5.OurresultsagreewillwiththereferencetimeTr=8318:63in[16].Wealsocomparetheripeningtimeusingseveralschemes,withlargertimestepsinTable3.6.Amongthethreemethods,BDF2isthemostt,andprovidesbetterestimatesofthetrueripeningtime,evenforlargert.Inthemostextremeinstanceoft=10,wenotethatthescheme(BE)predictsripeningtoosoon,andthatSDC2istoolate;butBDF2isstillgivesfairlyaccurateevolution.However,tocapturetheripeningmomentaccurately,westillrequiresmalltimestepst0:05),whichistooexpensiveforlongtimesimulations.Thus,weconsideradaptivetimesteppingforthesameproblem.InSection3.6,weexplainedthattheremightbeseveralwaystoapproximatethelocaltruncationerror(LTE).Withourschemes,weusethreetmethodsfortheLTEapproximation:RichardsonextrapolationbasedonBEscheme(inAlgorithm1);orBEcombinedwithBDF2(BE-BDF2)orcombinedwithSDC2(BE-SDC2).WenotethatwedonotconsiderSDIRK2,sinceithasalreadyprovedthatthisschemeistthantwoothersinSection3.3.4.WeimplementthreemethodswiththesamedpointresidualtoleranceNtol=1011,Nmaxit=600inAlgorithm1,butwithvariouserrortolerancetol.Theperformanceofthetime-adaptiveschemeisshowninTable3.7,whicharemoreaccurateandtthanprevioustime-steppingmethods.Inparticular,BE-BDF2combinationforLTEapproximationisthefastest,aswellasthemostaccurate,inthissimulation.Thephasefunctionu(x;t)obtainedbyouradaptivetimestepping(BE-BDF2)scheme83RipeningtimeBE8317:81BDF28318:70BDF38318:74SDC28318:99SDC38318:84Table3.5:Ripeningtimeof1DCHequationwithsmalltimestepsizet=0:01)andwiththeperiodicBC.timestepRipeningtimeTimes(s)108250:00313:38BE18296:00351:520:058311:851312:16109050:00717:26SDC218582:00692:650:058319:802779:61108290:00231:40BDF218303:00281:560:058315:751288:68Table3.6:Ripeningtimeof1DCHequationwithlargetimestepst=0:05;1;10)andwiththeperiodicBC.tolRipeningtimeTimes(s)1038292:54224:03BE1048276:08226:06(Richardson)1058308:23235:441038311:08233:38BE-SDC21048311:47239:581058312:91226:431038320:03184:71BE-BDF21048319:91196:821058317:87208:18Table3.7:Ripeningtimeof1DCHequationwithadaptivetimestepsizeandwiththeperiodicBC.84(a)u(x;0)(b)u(x;0:51)(c)u(x;3669:8)(d)u(x;7005:7)(e)u(x;8317:9)(f)u(x;8319:2)Figure3.5:Temporalevolutionofthe1DCHsolutionfromtheadaptivetimestep(BE-BDF2)scheme.isshowninFigure3.5.Theinitialstate3.5aquicklymovestothemetastablestate3.5b,andthenreachesthestablestate3.5fatwhichtwolayersaremergedtogetherafteraveryslowtimescale.ThissimulationalsoshowsthattheripeningeventhappensoveraveryfasttimescaleinFigure3.5e.InFigure3.6,wealsoplotthehistoryoftimestepsize,iterationcountateachtimelevel,anddiscreteenergyduringtimeevolutionsofournumericalsolutionwhichisobtainedbyadaptivetimesteppingmethod(BE-BDF2).AsshowninFigure3.6a,smalltimestepsareusedatearlystage(spinodalevolution)butincreasewhenthecoarseningprocessstarts,whichspeedupthesimulation.Also,ifiterationcountistoolarge(NitNmaxit1),werejectthesolutionandcomputeu(x;t)againwiththereducedtimestept.WeseethisbehaviorinFigure3.6b,wheretimestepsaredecreasedwhenverNitˇNmaxit.Wealsoobservein85(a)Timestepsize(b)Iterationcount(c)EnergyFigure3.6:Thetimestep,iterationcount,andenergyhistoryofouradaptivetimestep(BE-BDF2)schemefor1DCHequation.Figure3.6cthatadaptivetimesteppingdoesnotenergydecay.3.7.22DCahn-HilliardModelWenextsolvetheCHequation(3.1)intwospatialdimension.Withtheparameters=0:18;x=y=2ˇ128ˇ0:0491;Nmaxit=500;Ntol=107;weimplementtheinitialcondition(3.37),usingboththeandadaptive(BE-BDF2)timesteppingmethods.Weimplementtimesteppingbyusingtherder(BE),second-order(BDF2),andthethird-order(BDF3)methodswithasmalltimestept=0:01)andchecktheripeningtime.Thenumericalripeningtimeisatwhichu(ˇ2;ˇ2)changesfrompositivetonegative,becausethepoint(ˇ2;ˇ2)isthecenterofthesmallercircularregionofu=1,whichwillbeconsumedbythelargeroneafteralongtimescale.BasedonnumericalresultsinTable3.8,wecansuchripeningtimeasTr=80:07.Eachtimesteppingmethodcombinedwithdimensional-splittingguaranteesthatthe86TimesteppingRipeningtimeTotaliterationTimes(s)BE80:04167,0862025:96t=0:01BDF280:07149,5001883:72BDF380:07142,2401793:91Table3.8:Numericalripeningtimeofthe2DCHequationwithxedtimesteppingmethodst=0:01,Ntol=107)andwiththeperiodicBCs.(a)Energy(BE)(b)Energy(BDF2)(c)Energy(BDF3)Figure3.7:Thediscreteenergyofedtimesteppingschemes((a)BE,(b)BDF2,and(c)BDF3)for2DCHequation.Thecommonparametersareused:t=5,andx=y=0:05solutionisdecayinginenergy,whichshowninFigure3.7,evenwheretheCFLis100,t=5,andx=y=0:05).Bytheseresults,westronglybelieveourhigher-orderdimensionalsplittingMOLTmethodsareunconditionalgradientstableinpractice.Inaddition,itisimportanttonotethatthatraisingtheorderofthescheme(temporalaccuracy)reducesthenumberofiterationpertimestep,especiallyattheripeningmoment,asshowninFigure3.8.Thus,theoverallcomputationtimeofhigher-ordermethodislowerthanlow-orderschemes,summarizedinTable3.8.Wealsoimplementthesametimesteppingmethodswithalargertimestept=0:1.AsshowninFigure3.8b,eachofiterationcountspertimestepishigherthanthecorrespondingonewithsmallertimestep(inFigure3.8a).Webelievethatthisisreasonableresultsbecauseoftheaccuracyrequirement.Moreover,thehigher-ordertimesteppingschemeexhibitsthebetterlikethesmallertime87(a)t=0:01(b)t=0:1Figure3.8:Iterationcountof2DCHsolutionobtainedbythetimestepppingmethod(a)t=0:01),and(b)t=0:1).BlacklineisobtainedbyBackwardEuler(BE)scheme;blueisBDF2;andtheredoneisBDF3atbothplots.stepcase.Withthesetests,westronglybelievethathigher-ordertimesteppingmethodispreferabletoresolvetheinterfacialstructureofthesetypesoftheproblems,regardlessoftimestepsizes.However,t=0:1orevent=0:01isstilllargeforthisCHproblem,especiallyduringspinodalcomposition,sothatitispossibletodampouthighfrequencycontributionsintheevolvingwhichissomewhatunphysical.Thus,wenextsimulatethesameproblemwithadaptivetimesteppingmethod(BE-BDF2).Wecomputethenumericalripeningtime,totaliterationnumberduringtemporalevolution,andthecomputationaltimeofouradaptivetimesteppingmethod,showninTable3.9.Whenwecomparethisresultwiththetimesteppingmethods,theadaptiveschemeisbetterwithrespecttotheaccracy,asexpected.Thereasonwhythecomputationaltimeofadaptiveschemeisalittlelargerthanhigher-orderschemesisitusesverysmalltimestepsizeduringspinodalevolution.Thus,werealizethattheadaptivetimesteppingmethodisnecessaryfortheCHmodelinmultiple88TimesteppingRipeningtimeTotaliterationTimes(s)BE79:8084,5221125:69t=0:1BDF280:1068,322844:39BDF380:1061,548770:67AdaptivetimeBE-BDF280:08366,085918:53Table3.9:Numericalripeningtimeofthe2DCHequation(theperiodicBCs):Comparisonbetweentheedtimesteppingmethodst=0:1,Ntol=107)andtheadaptivetimesteppingmethod(BE-BDF2).spatialdimension,aswell.Like1Dcase,thetimestepsize,iterationcount,andenergyhistoryarepresentedinFigure3.9.OurresultsagreewiththoseofSection3.7.1,inthatlargertimestepsizesareusedwherethetransitionlayersarevaryingslowly,andsmallerstepsareusedwherethelayersvaryrapidly(seeFigure3.9a),sothatwebelieveourschemecapturesvarioustimescalesof2DCHsolutionaccurately.TheenergyinFigure3.9cindicatesthattherearetwosharptransitionsintheenergyE;earlyon,andatripening.Finally,wealsoobservethedesiredenergydecayingpropertyofournumericalsolution.Thecontourplotsoftimeevolutionofthephasefunctionu(x;y;t)areshowninFigure3.10.Theinitialstates3.10aquicklyreachesthemetastablestate,whereweseetwocircularformations.Eventuallythelargeroneabsorbsthesmaller,althoughoveraverylongtimescale,shownin3.10c-3.10e.Thestateisshownin3.10f,wherethelargerregionhasfullyconsumedthesmaller.Inallplots,thetotalvolumelookstobepreserved.ThecomputationalripeningtimeisTrˇ80:0834inthissimulation.3.7.32DvectorCahn-HilliardModelWenowapplyourandsecond-ordertimesteppingmethodsfromSection3.5tothe2DvectorCHsystem(3.3)combinedwithouradaptivetimestepping(BE-BDF2)strategy.In89(a)Timestepsize(b)Iterationcount(c)EnergyFigure3.9:Thetimestep,iterationcount,andenergyhistoryofouradaptivetimestep(BE-BDF2)schemefor2DCHequation.(a)u(x;y;0)(b)u(x;y;0:5083)(c)u(x;y;15:807)(d)u(x;y;70:015)(e)u(x;y;80:083)(f)u(x;y;90)Figure3.10:Temporalevolutionofthe2DCHsolutionwiththeinitial(3.37)..90thisSection,weobservethelongtimebehaviorofthephasefunctionu=(u1;u2),usingthesameinitialconditionforu1andu2inSection3.5.1.Withtheparameters=0:32;x=y=2ˇ64ˇ0:0982;Nmaxit=400;Ntol=106;tol=104weimplementtheadaptivetimesteppingmethod(BE-BDF2)forvectorCH(VCH)model.Firstofall,thetimestep,numberofnonlineariterationandenergyhistoryresultsof2DVCHsolutionarepresentedinFigure3.11.OurresultsagreewithpreviousSections,wealsoobservethedesiredenergydecayinFigure3.11c.Inaddition,thecontourplotsofcos(arg(u1+iu2))areshowninFigure3.12.Insteadofplottingu1andu2separately,[16]wetheanglearg(u1+iu2)atatriplejuction,andplotcos()toavoidanydiscontinuities.Wefollowthisbenchmarkplotusingournumericalsolution.AftersomeinitialripeninginFigure3.12b,theinterfacesdividingthethreestatesu=zi(i=1;2;3)formaroundT=0:5.InFigure3.12c-3.12f,twoofthevalueshavecos(2ˇ=3)(lightblueintheplots)butseparatedbydarkbluelines.Thentheripeningprocessbegins,andoccursoveralongtimescale,endingaroundT=24:286inFigure3.12e.Again,webelievethatthevolumeispreservedoveralltimesteps.Foramoreinvolveddiscussionofthissimulation,werefertheinterestedreaderto[16].Ourgoalhereistoreproducethesameresultsinthereference,butbasedonourdimensional-splitMOLTschemeandthatourschemecancapturethecorrectripeningbehaviorinmultiplespatialdomain.Moreover,thenumericalresultswhatwehavedoneinthisSectionprovethatourschemecanbeextensibletovectorsystem.Futureworkwillinvestigatetheextensionofthesesecond-andthird-ordertimesteppingschemesofvectormodeltoapplythemulticomponentbilayer91structuresariseinmembraneinbiology[51].3.7.42DsixthorderModelInthisSection,weintrouduceoneofourongoingresearchprojects.Wealsoconsiderthesixthordereldproblem[16]:ut=h(2f0(u)+2)(2uf(u))i;f(u)=u3u:(3.46)whereandaregivenpositiveconstants.(Thesignofisimportantsinceminuscanformthephaseinterface)ThisproblemismotivatedbythefunctionalizedCahn-Hilliard(FCH)equation[29,21,37,40]whichmodelsinterfacialenergyinamphiphilicphase-separatedmixtures.IntheoriginalFCHmodel,1>0and22R,whicharefunctionalizationtermthatareanalogoustothesurfaceandvolumeenergiestypicalofmodelsofchargedsolutes,butwehavetothecase=1=2.Later,wewillconsiderthesolutionoforiginalFCHequationwithourMOLTformulation.WiththesimilarapproachestoCHtypeequations,weintroducethetransformedvariablev=u+1andsubstituesin(3.46),vt=h(2f0(v)+2)(2vf(v))i;f(v)=v33v2+2v:(3.47)ApplythebackwardEuler(BE)schemefortimediscretizationof(3.47),I43vn+1=vnt2fn+12(f0v)n+1+(ff0)n+12t2vn+1fn+1:92(a)Timestepsize(b)Iterationcount(c)EnergyFigure3.11:Thetimestep,iterationcount,andenergyhistoryofVCHsolution.(a)t=0(b)t=0:02503(c)t=0:5083(d)t=5:0848(e)t=24:286(f)t=30:021Figure3.12:Temporalevolutionofthe2DCHvectorsolutionwiththeinitial(3.37)foru1and(3.44)foru2.Contoursofcos(argu1+iu2)areplotted.Twophasesu=z2andu=z3havesamecosinevalues,12(lightblueintheplots)butareseparatedbydarkbluelines.93Hence,wenowshouldinvertthe6thorderoperatortosolvevn+1.OnecaninvertthisbycompletingthecubesuchthatI3p43vn+1=vnt2fn+12(f0v)n+1+(ff0)n+12t2vn+1fn+133p43(3p4)22vn+1:(3.48)Wecannowsolvevn+1byapplyingthetripleinversionofourmoHelmhotlzoperatorL=I3p4bythesameprocedureinSection3.4.1.Forsimplicity,weapplytheBEscheme(3.48)withthetimestept=0:1tosolvethe6thorderproblem(3.47)(=1).Bystartingwiththesameinitialcondition(3.37)(=0:18),andwiththefollowingparameters:x=y=2ˇ128ˇ0:05;Nmaxit=200;Ntol=106:Thecontourplotsofthetemporalevolutionsofournumericalsolutionu(x;y;t)of(3.46)areshowninFigure3.14.Aftersolutionvariesrapidly,thensolutionseemstomoveveryslowly.Asexpected[16],webelievethatournumericalsolutionreachestothesteadystatearoundtˇ320,whichhasformedtheregulararrayinFigure3.14fandthecorre-spondingdiscreteenergydecreasedatthatmomentinFigure3.13b.Futureworkwillconsidertheswitching-schemeforthesixth-orderproblem,andextendtotheadaptivetimesteppingstrategy,similarwiththepreviouslypresentedCHsolutions.94(a)Iterationcount(b)EnergyFigure3.13:Theiterationcountandenergyhistoryofthesolutionofsixth-orderequation(3.47).(a)u(x;y;0)(b)u(x;y;2)(c)u(x;y;40)(d)u(x;y;150)(e)u(x;y;300)(f)u(x;y;500)Figure3.14:Temporalevolutionofthe2Dsixthordermodel'ssolutionwiththeinitial(3.37).95Chapter4ConclusionandFutureworkInthiswork,wehaveproposedtheMethodoflinestranspose(MOLT)formulationfornonlinearproblemssuchasAllen-CahnandCahn-Hilliardequations.Thereasonwehavedoneisthatthesemi-discreteformulationgivesbetterfavorablestabilitypropertythanatraditionalMethodoflines(MOL)formulation.Then,wehaveprovedthesuccessiveconvolution(resolventexpansion)couldachievearbitraryorderintime,givingusgreattimeaccuracy.Moreover,wehavecomparedmoretraditionaltimesteppingalgorithms,suchasbackwardformula(BDF),singlydiagonallyimplicitRungekutta(SDIRK),andspectraldeferredcorrection(SDC),toachievehighorderofaccuracyforCahn-Hilliardmodels.Furthermore,theadaptivetimesteppingisstillfavorabletocapturethescaleseparationofCHmodel.Sp,smalltimestepsareusedatthespinodalphase,butlargestepsareusedatthecoarseningprocessinwhichphase-separateddomainsmergeintolargerdomains.Inordertohandlethespatialcomponentoftheproblems,typicallythefastFouriermethodarechosen,butthisonlyguaranteestheperiodicboundaryconditions.Instead,wehaveintroducedadimensional-splitkernel,whichistO(N),andmatrix-freescheme.Notonlywecanaddressotherboundaryconditions,wecanraisearbitraryorderofaccuracyinspace.Thisalsopromisestheparallelncyformodernparallelmulticorecomputingbecauseofthedimensionalsplittingstrategy.Especially,thedimensionalsplittingalsohasaframeworkthatwecanextendour1D96ToolsResultsMOLTdiscretizationStabilityinsemi-discretesettingTimeSuccessiveconvolutionHigh-orderintimesteppingBDF/SDIRK/SDCComparestrategiesTime-adaptiveCapturevarioustimescalesFlexibilitywithBCsSpatialFastArbitraryorderofaccuracyinspaceschemeConvolutionFastO(N)computationParallelMulti-DDimensionalEasyextensiontomulti-DschemeschemesplittingHigh-orderspaceandtimeTable4.1:Summarizationofnumericalresultsofpresentedscheme,dimensionalsplitMOLTscheme.schemetomultiplespatialdimensions.In2D,wehaveachievedthesameresultswith1D,i.e.t,canachievehigh-orderofaccuracy,andaddressvariousboundaryconditions.Onemightwonderthesplittingerrorintroducedbythedimensionalsplitting,however,wehavealsoaddressedthesplittingerrorinouriterativesolver.WesummarizethenumericalstrategieswehavedoneinthisworkandresultingpositivewhatwehavegotintheTable4.1,whichiseasytoseeaaglanceourworks.Infuturework,parallelimplementationwillbethefocusofourupcomingworkinthisarea.Moreover,wewillinvestigateemployingourschemetothehigher-orderderivativemodels,suchas(2Dand3D)FCHmodelsin[40].Inthiswork,wehavedonetimesteppingmethodforthesixthorderparabolicequations,andthenextgoalshouldemployadaptivetimestepping,whichcanallowincreasedtimestepstosolvesuchhigher-orderderivativemodels.Wealsowanttomakeuseoftheimplicitsolversabilitytohandlecomplexboundarygeometriesofvariousmodels.WehavenotyettriedcomparisonstudybetweenotherboundaryintegralapproximationtosolvethemoHelmholtzequationsintheMOLTsemi-discreteschemes,suchasFFM97inSection1.2.2.1orTreecodealgorithm[45].Besidesourdimensionalsplittingalgorithm,wewillimplementotherboundaryintegralsolvers,andcomparewhichsolversgainciencyandaccuracyinparallelcomputing.98BIBLIOGRAPHY99BIBLIOGRAPHY[1]L.AbadiasandP.J.Miana.c0-semigroupsandresolventoperatorsapproximatedbylaguerreexpansions.arXivpreprintarXiv:1311.7542,2013.[2]R.Alexander.Diagonallyimplicitrunge-kuttamethodsforode's.SIAMJournalonNumericalAnalysis,14(6):1006{1021,1977.[3]S.M.AllenandJ.W.Cahn.Amicroscopictheoryforantiphaseboundarymotionanditsapplicationtoantiphasedomaincoarsening.ActaMetallurgica,27(6):1085{1095,1979.[4]D.M.Anderson,G.B.McFadden,andA.A.Wheeler.terfacemethodsinmechanics.Annualreviewofmechanics,30(1):139{165,1998.[5]L.BronsardandF.Reitich.Onthree-phaseboundarymotionandthesingularlimitofavector-valuedginzburg-landauequation.ArchiveforRationalMechanicsandAnalysis,124(4):355{379,1993.[6]O.P.BrunoandM.Lyon.High-orderunconditionallystablefc-adsolversforgeneralsmoothdomainsi.basicelements.JournalofComputationalPhysics,229(6):2009{2033,2010.[7]J.W.Cahn.Onspinodaldecomposition.Actametallurgica,9(9):795{801,1961.[8]J.W.CahnandJ.E.Hilliard.Freeenergyofanonuniformsystem.i.interfacialfreeenergy.TheJournalofchemicalphysics,28(2):258{267,1958.[9]J.CarrandR.L.Pego.Metastablepatternsinsolutionsofut=2uxx-f(u).Commu-nicationsonpureandappliedmathematics,42(5):523{576,1989.[10]M.Causley,A.Christlieb,Y.Gu˘clu,andE.Wolf.Methodoflinestranspose:anta-stablesolverforwavepropagation.arXivpreprintarXiv:1511.01013,2015.[11]M.Causley,A.Christlieb,B.Ong,andL.VanGroningen.Methodoflinestranspose:Animplicitsolutiontothewaveequation.MathematicsofComputation,83(290):2763{2786,2014.100[12]M.F.Causley,H.Cho,A.J.Christlieb,andD.C.Seal.Methodoflinestranspose:Highorderl-stableo(n)schemesforparabolicequationsusingsuccessiveconvolution.SIAMJournalonNumericalAnalysis,54(3):1635{1652,2016.[13]M.F.CausleyandA.J.Christlieb.Higherordera-stableschemesforthewaveequa-tionusingasuccessiveconvolutionapproach.SIAMJournalonNumericalAnalysis,52(1):220{235,2014.[14]L.ChenandJ.Shen.Applicationsofsemi-implicitfourier-spectralmethodtophaseequations.ComputerPhysicsCommunications,108(2):147{158,1998.[15]L.-Q.Chen.Phasemodelsformicrostructureevolution.Annualreviewofmaterialsresearch,32(113{140),2002.[16]A.Christlieb,J.Jones,K.Promislow,B.Wetton,andM.Willoughby.Highaccuracyso-lutionstoenergygradientwsfrommaterialsciencemodels.JournalofComputationalPhysics,257:193{215,2014.[17]A.Christlieb,C.Macdonald,B.Ong,andR.Spiteri.Revisionistintegraldeferredcorrectionwithadaptivestep-sizecontrol.CommunicationsinAppliedMathematicsandComputationalScience,10(1):1{25,2015.[18]A.Christlieb,K.Promislow,andZ.Xu.Ontheunconditionallygradientstableschemeforthecahn-hilliardequationanditsimplementationwithfouriermethod.Commun.Math.Sci,11:345{360,2013.[19]A.J.Christlieb,Y.Liu,andZ.Xu.Highorderoperatorsplittingmethodsbasedonanintegraldeferredcorrectionframework.JournalofComputationalPhysics,294:224{242,2015.[20]S.M.CoxandP.C.Matthews.Exponentialtimeforsystems.JournalofComputationalPhysics,176(2):430{455,2002.[21]S.DaiandK.Promislow.Geometricevolutionofbilayersunderthefunctionalizedcahn{hilliardequation.InProceedingsoftheRoyalSocietyofLondonA:Mathematical,PhysicalandEngineeringSciences,volume469,page20120505.TheRoyalSociety,2013.[22]J.Douglas,Jr.Onthenumericalintegrationof@2u@x2+@2u@y2=@u@tbyimplicitmethods.Journalofthesocietyforindustrialandappliedmathematics,3(1):42{65,1955.101[23]A.Dutt,L.Greengard,andV.Rokhlin.Spectraldeferredcorrectionmethodsforordi-narytialequations.BITNumericalMathematics,40(2):241{266,2000.[24]K.-J.EngelandR.Nagel.One-parametersemigroupsforlinearevolutionequations,volume194.SpringerScience&BusinessMedia,2000.[25]L.C.Evans.PartialentialEquations:SecondEdition.AmericanMathematicalSociety,2010.[26]D.J.Eyre.Anunconditionallystableone-stepschemeforgradientsystems.Unpublishedarticle,1998.[27]G.FairweatherandA.Mitchell.Anewcomputationalprocedureforadimethods.SIAMJournalonNumericalAnalysis,4(2):163{170,1967.[28]P.C.Fife.Mathematicalaspectsofreactingandsystems,volume28.SpringerScience&BusinessMedia,2013.[29]N.Gavish,J.Jones,Z.Xu,A.Christlieb,andK.Promislow.Variationalmodelsofnetworkformationandiontransport:applicationstopionomermem-branes.Polymers,4(1):630{655,2012.[30]J.W.Gibbs,H.A.Bumstead,andW.R.Longley.ThecollectedworksofJ.WillardGibbs,volume1.Longmans,GreenandCompany,1928.[31]L.GreengardandV.Rokhlin.Afastalgorithmforparticlesimulations.Journalofcomputationalphysics,73(2):325{348,1987.[32]L.F.GreengardandJ.Huang.Anewversionofthefastmultipolemethodforscreenedcoulombinteractionsinthreedimensions.JournalofComputationalPhysics,180(2):642{658,2002.[33]V.GrimmandM.Gugat.Approximationofsemigroupsandrelatedoperatorfunctionsbyresolventseries.SIAMJournalonNumericalAnalysis,48(5):1826{1845,2010.[34]A.Iserles.Acourseinthenumericalanalysisofentialequations.Number44.Cambridgeuniversitypress,2009.[35]J.JiaandJ.Huang.Krylovdeferredcorrectionacceleratedmethodoflinestransposeforparabolicproblems.JournalofComputationalPhysics,227(3):1739{1753,2008.102[36]S.Jiang,L.Greengard,andS.Wang.cientsum-of-exponentialsapproximationsfortheheatkernelandtheirapplications.AdvancesinComputationalMathematics,41(3):529{551,2015.[37]J.S.Jones.Developmentofafastandaccuratetimesteppingschemeforthefunction-alizedCahn-Hilliardequationandapplicationtoagraphicsprocessingunit.PhDthesis,MICHIGANSTATEUNIVERSITY,2013.[38]J.Kacur.MethodofRotheinevolutionequations.Springer,1986.[39]J.P.KeenerandJ.Sneyd.Mathematicalphysiology,volume1.Springer,1998.[40]N.KraitzmanandK.Promislow.Anoverviewofnetworkbifurcationsinthefunction-alizedcahn-hilliardfreeenergy.InMathematicsofEnergyandClimateChange,pages191{214.Springer,2015.[41]V.KrinskyandA.Pumir.Modelsofofcardiactissue.Chaos:AnInter-disciplinaryJournalofNonlinearScience,8(1):188{203,1998.[42]M.C.A.KropinskiandB.D.Quaife.Fastintegralequationmethodsforthemohelmholtzequation.JournalofComputationalPhysics,230(2):425{434,2011.[43]L.D.LandauandE.Lifshitz.Onthetheoryofthedispersionofmagneticpermeabilityinferromagneticbodies.Phys.Z.Sowjetunion,8(153):101{114,1935.[44]H.G.LeeandJ.-Y.Lee.Asemi-analyticalfourierspectralmethodfortheallen{cahnequation.Computers&MathematicswithApplications,68(3):174{184,2014.[45]K.LindsayandR.Krasny.Aparticlemethodandadaptivetreecodeforvortexsheetmotioninthree-dimensionalw.JournalofComputationalPhysics,172(2):879{907,2001.[46]F.LiuandJ.Shen.Stabilizedsemi-implicitspectraldeferredcorrectionmethodsforallen-cahnandcahn-hilliardequations.Math.MethodsAppl.Sci,2013.[47]C.LubichandR.Schneider.Timediscretizationofparabolicboundaryintegralequa-tions.NumerischeMathematik,63(1):455{481,1992.[48]M.LyonandO.P.Bruno.High-orderunconditionallystablefc-adsolversforgeneralsmoothdomainsii.elliptic,parabolicandhyperbolicpdes;theoreticalconsiderations.JournalofComputationalPhysics,229(9):3358{3381,2010.103[49]S.P.N˝rsett.One-stepmethodsofhermitetypefornumericalintegrationofsystems.BITNumericalMathematics,14(1):63{77,1974.[50]D.OlmosandB.D.Shizgal.Pseudospectralmethodofsolutionoftheugh{nagumoequation.MathematicsandComputersinSimulation,79(7):2258{2278,2009.[51]K.PromislowandQ.Wu.Geometricevolutionofquasi-bilayersinmulticomponentfunctionalizedcahn-hilliardequation.arXivpreprintarXiv:1510.08467,2015.[52]W.E.Schiesser.Thenumericalmethodoflines:integrationofpartialentialequa-tions.Elsevier,2012.[53]D.C.Seal,Y.Gu˘clu,andA.J.Christlieb.High-ordermultiderivativetimeintegratorsforhyperbolicconservationlaws.JournalofComputing,60(1):101{140,2014.[54]J.ShenandX.Yang.Numericalapproximationsofallen-cahnandcahn-hilliardequa-tions.DiscreteContin.Dyn.Syst,28(4):1669{1691,2010.[55]J.M.StarobinandC.Starmer.Commonmechanismlinksspiralwavemeanderingandwave-front{obstacleseparation.PhysicalReviewE,55(1):1193,1997.[56]G.TierraandF.Guillalez.Numericalmethodsforsolvingthecahn{hilliardequationanditsapplicabilitytorelatedenergy-basedmodels.ArchivesofComputationalMethodsinEngineering,2(269{289):2015,22.[57]M.R.Willoughby.High-ordertime-adaptivenumericalmethodsfortheallen-cahnandcahn-hilliardequations.Master'sthesis,UniversityofBritishColumbia,2011.[58]E.M.Wolf.Aparticle-in-cellmethodforthesimulationofplasmasbasedonanuncon-ditionallystablewaveequationsolver.PhDthesis,MICHIGANSTATEUNIVERSITY,2015.[59]F.Xie,Z.Qu,J.N.Weiss,andA.el.Coexistenceofmultiplespiralwaveswithindependentfrequenciesinaheterogeneousexcitablemedium.PhysicalReviewE,63(3):031905,2001.104