lecture13[1]

42NUMERICALSOLUTIONOFPDES6.Efficientsolutionofthelinearsystemsarisingfromfiniteelementdiscretization6.1.Optimizationmethods.WehaveshownthatthefiniteelementdiscretizationofPois-son’sequationleadstothesolutionofalinearsystemAx=b,inwhichAisasymmetricmatrix.ItisalsoeasytocheckthatAisapositivedefinitematrix,i.e.,xTAx0forx6=0.Forsuchaproblem,thesolutionofthissystemisalsotheminimizerofthefunctionalφ(x)=12xTAx−xTb.Notetheminimumwilloccurwhere∇φ(x)=0.But∇φ(x)=Ax−b,sothesolutionoftheminimizationproblemisthesolutionofthelinearsystemofequations.Atypicalminimizationalgorithmistolet{pk}k≥0beasetofsearchdirectionsand{αk}k≥0asetofscalarsanddefineaniterationxk+1=xk+αkpk.Thesimplestexampleisthemethodofsteepestdescent,inwhichwechoosepk=−∇φ(xk)=−[Axk−b].Todeterminethebestchoiceofαk,wethenminimizeφ(xk+αkpk)withrespecttoαk,consideringxkandpknowfixed.Sinceφ(xk+αkpk)=12(xk)TAxk+2αk(pk)TAxk+α2k(pk)TApk−xTb−αkpTb,minimizingwithrespecttoαkgives:(pk)TAxk+αk(pk)TApk−(pk)Tb=0,i.e,αk=(pk)T(b−Axk)(pk)TApk=(pk)Tpk(pk)TApk.Ifwemakethesimplerchoice,αk=αforallk,thenwegettheiterationxk+1=xk−α[Axk−b]=[I−αA]xk+αb.IfweletxdenotetheexactsolutionofAx=b,thenwegettheerrorequationx−xk+1=x−[I−αA]xk−αb=[I−αA](x−xk)+αAx−αb=[I−αA](x−xk).Iteratingthisequation,wefindthatx−xk=[I−αA]k(x−x0).Awellknownresultfromlinearalgebrasaysthatthisiterationwillconvergeforallx0∈Rnifandonlyifρ(I−αA)1,whereifMisann×nmatrixwitheigenvaluesμi,thenρ(M)=maxi|μi|.NowifλisaneigenvalueofA,then1−αλisaneigenvalueofI−αA(withthesameeigenvector).Hence,forconvergence,weneed−11−αλ1foralleigenvaluesλofthematrixA.SinceAispositivedefinite,allitseigenvaluesarepositive,sowerequire0α2/λ,i.e.,0α2/ρ(A).Todeterminetheoptimalchoiceoftheparameterα,weminimizethenormoftheiterationmatrixI−αA.IfweconsiderkI−αAk2,thensinceAisassumedsymmetric,soisI−αA.NUMERICALSOLUTIONOFPDES43Hence,kI−αAk2=ρ(I−αA)=maxi|1−αλi|,whereλiaretheeigenvaluesofA.SinceAispositivedefinite,wehavethat0λ1≤...≤λn.Thenmaxi|1−αλi|=max{|1−αλ1|,|1−αλn|}andthismaximumwilloccurwherethesetwoquantitiesareequal,i.e.,1−αλ1=αλn−1.Hence,theoptimalvalueisα=2/(λ1+λn).Inthiscase,ρ(I−αA)=1−2λ1λ1+λn=λn−λ1λn+λ1=(λn/λ1)−1(λn/λ1)+1.Letκ=kAk2kA−1k2betheconditionnumbermeasuredinthek·k2norm.SinceAissymmetricandpositivedefinite,kAk2=ρ(A)=λn.SincetheeigenvaluesofA−1arethereciprocalsoftheeigenvaluesofA,kA−1k2=ρ(A−1)=1/λ1.Hence,κ=λn/λ1.Thus,ρ(I−αA)=(κ−1)/(κ+1),andwehaveprovedthefollowingresult.Theorem9.IfAissymmetricandpositivedefinite,thentheiterationschemedefinedbyxk+1=[I−αA]xk+αb,withα=2/(λ1+λn)satisfies:kx−xkk2≤κ−1κ+1kkx−x0k2,κ=λmax(A)/λmin(A).ForthesolutionofPoisson’sproblembystandardfiniteelements,wecanshowthatthereisaconstantindependentofhsuchthatκ(A)≤c2h−2.Thus,implementingthisiterationinitspresentformleadstoasmallreductioninerror(1−O(h2))andslowconvergence.Togetamorepreciseunderstandingofwhatthemethodisdoing,weconsideraneigen-functionexpansionoftheerror,i.e.,wesupposethatAφi=λiφi,where{φi}Ni=1areasetoforthonormaleigenvectorsofA.Wethensetek=x−xkandwritee0=NXi=1[(e0)Tφi]φi.Supposewechooseα=1/λN,thelargesteigenvalueofA.Thenek=[I−αA]ke0=NXi=1[(e0)Tφi](1−λi/λN)kφi.Nowforlargeeigenvalues1−λi/λNissmall,sothehighfrequencycomponentsoftheerroraredampedoutquickly,whileforsmalleigenvalues1−λi/λN≈1,andthereisnotmuchdecayintheerrorandsothelowfrequencycomponentsarenotchangedmuch.Thus,afewiterationsofthismethodhastheeffectof“smoothing”theerror.Weshallcomebacktothisideainalaterlecture.6.2.Conjugate-Gradientmethod(CG).Abetterchoiceofsearchdirections{pk}istochoosethemtobeA-orthogonal,i.e,tosatisfy(pj)TApi=0fori6=j.Inthiscase,thebestchoiceoftheαkaregivenbyαk=(pk)T[b−Axk](pk)TApk=(pk)T[b−Ax0](pk)TApk.44NUMERICALSOLUTIONOFPDESTheCGmethodgeneratesthedirectionspkrecursivelyusingtheGram-Schmidtorthogo-nalizationprocess(seereferencesforprecisedescription).FortheCGmethod,wegetthefollowingerrorestimate:kx−xkkA≤2√κ−1√κ+1kkx−x0kA,wherekxk2A=xTAx.Sincenow√κenters,thereductionislike1−O(h),betterthanbefore,butstillslow.Inpractice,oneusestheideaofpreconditioning.InsteadofsolvingthesystemAx=b,wesolvethesystemBAx=Bb,whereBisanapproximationtoA−1thatiseasilycomputable.ThentherateofconvergencedependsontheconditionnumberofBAinsteadofA.IfBisagoodapproximationtoA−1,thenBA≈I,andsoκ(BA)willbecloseto1,andwewillgetasubstantialerrorreductionateachiteration.