Table Of Content

N M I UMERISCHE ATHEMATIK Wintersemester 2009/2010 Leibniz Universität Hannover PROF. DR. MARC STEINBACH SkriptzurVorlesung Prof.Dr.MarcSteinbach InstitutfürAngewandteMathematik LeibnizUniversitätHannover Welfengarten1 30167Hannover [email protected] LiteraturempfehlungenzurVorlesung P.Deuflhard,A.Hohmann. NumerischeMathematikI. DeGruyter,Berlin,1993. M. Hanke-Bourgeois. Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens. Teubner,Stuttgart,2002. A.Quarteroni,R.Sacco,F.Saleri. NumerischeMathematik1. Springer,Berlin,2002. A.Quarteroni,R.Sacco,F.Saleri. NumerischeMathematik2. Springer,Berlin,2002. Kapitel 0 Aufgaben und Ziele der Numerischen Mathematik Modellgestützte Analyse, Simulation und Optimierung von Vorgängen in Naturwissenschaften, Technik, Wirtschaft usw. (Physik, Chemie, Biologie; Baumechanik, Strömungsmechanik, Elektrotechnik, Mecha- tronik;Finanzwesen,Aktien-,Waren-undEnergiemärkte): 1. MathematischeModellbildungmittelsexperimentellverifizierterGesetzeführtaufeinemathema- tische Problemstellung, z.B. eine Minimierungsaufgabe. In anspruchsvollen technischen Anwen- dungentritthierbeiofteinAnfangs-oderRandwertproblemfüreineDifferentialgleichungauf.Ein brauchbares mathematisches Modell sollte eine eindeutige Lösung besitzen und diese sollte stetig vondenvorgegebenenDaten(z.B.denAnfangs-oderRandwerten)abhängen. 2. Die Lösung solcher Aufgaben lässt sich in der Regel nicht geschlossen angeben, auch nicht unter Einsatz vonComputeralgebra-Software. Stattdessenmüssennumerische Verfahren verwendet werden, mit denen eine Näherungslösung berechnet wird. Aufgabe der Numerischen Mathematik istdieEntwicklungundAnalysenumerischerVerfahren.MitnumerischenRechenmethodenwerden aufdemComputerapproximativeZahlenwertefürdieLösungmathematischerAufgabenstellungen bestimmt. Diese unterscheiden sich grundlegend von symbolischen Rechenverfahren der Compu- teralgebra, bei denen mathematische Regeln (Umformungen, Differentiations- und Integrationsre- gelnusw.)zurManipulationvonFormelausdrückenverwendetwerden. 3. Visualisierung, d.h. geeignete graphische Darstellung der erhaltenen Näherungslösung, um eine VorstellungvondensimuliertenVorgängenzugewinnenundevtl.dasModellunddienumerischen Verfahrenzuverbessern. Als Wegweiserdurch die in dieser Vorlesung behandelten numerischen Verfahren wird uns das folgende BeispieleinesmathematischenModellsdienen.MitdiesemModelllassensichdieverschiedenenGrund- aufgabendernumerischenMathematikimLaufederVorlesungverdeutlichen. Beispiel(ElastischeBiegungeinesBalkens). Ziel ist die Berechnung der Durchbiegung eines an beiden Enden eingespannten Balkens aufgrund einer aufliegendenLast. y 6 f(x) L x - ??????????? Gesuchtistdie Funktion y(x)für 0 6 x 6 L, welchedie Auslenkungdes Balkensbeschreibt. Die Aus- lenkungistphysikalischdadurchbestimmt,dasssichderBalkenineinemGleichgewichtszustandbefindet, 2 NumerischeMathematikI in dem die potentielle Energie einen minimalen Wert annimmt. Für die potentielle Energie des Balkens lässtsich,fürdenFall„kleiner“Auslenkungeny(x),ausdenGesetzenderKontinuumsmechanikfolgende Formelherleiten: EJ L L W(y)= y′′(x)2dx+ f(x)y(x)dx. 2 Z0 Z0 innere(Biege-)Energie Lastpotential Dabei ist E der sogenannte Elastizität|smodul,{zeine vom} M|aterial{azbhängi}ge positive Zahl. Ist s die Di- cke des Balkens, dann bezeichnet J = s4/12 das zugehörige Biegemoment. Die Funktion f schließlich bezeichnetdieaufdenBalkenwirkendeLastverteilung. Die Tatsache, dass der Balken an den beiden Enden eingespannt ist, wird durch die Randbedingungen y(0) = 0, y′(0) = 0(amlinken Ende)sowiey(L) = y′(L) = 0(amrechtenEnde) berücksichtigt.Das mathematischeModellfürdieBestimmungdesGleichgewichtszustandeseinesbelastetenBalkensbesteht nuninderfolgendenMinimierungsaufgabe: BestimmeunterallenzweimalstetigdifferenzierbarenFunktioneny C2[0,L],welchedieRandbedingun- ∈ geny(0)=y′(0)=0undy(L)=y′(L)=0erfüllen,diejenige,diedasfolgendeFunktionalminimiert: EJ L L W(y)= y′′(x)2dx+ f(x)y(x)dx. 2 Z0 Z0 Die Lösung dieses Minimierungsproblems lässt sich allerdings nur in seltenen Spezialfällen geschlossen angeben. Stattdessen ist man darauf angewiesen, unter Einsatz des Computers eine hinreichend genaue ApproximationmitnumerischenMethodenzuberechnen. NatürlichkannmandazunichtmitdemunendlichdimensionalenVektorraumderzweimalstetigenFunk- tionen, deren Funktionswert und Ableitung am linken und rechten Intervallende verschwinden, arbeiten. Stattdessen müssen diese Funktionen durch einfachere angenähert werden, die sich mit endlich vielen Werten beschreiben lassen, beispielsweise durch Polynome von vorgegebenem Grad. Dies führt uns auf dieThematikderInterpolationvonFunktionen,dieimfolgendenKapitelbehandeltwird.FürdieInte- grale,diebeiderAuswertungdesEnergiefunktionalsW(y)zuberechnensind,lässtsichinderRegelkeine expliziteFormelfürdieStammfunktionangeben.StattdessenmusshierzuaufMethodendernumerischen Integrationzurückgegriffenwerden,dieGegenstanddeszweitenKapitelsseinwerden. Danach folgt ein Kapitel mit grundsätzlichen Untersuchungen zu der Frage, inwieweit man sich auf die mittels numerischen Methoden berechneten Näherungslösungen überhaupt verlassen kann. Der Begriff der Kondition eines mathematischen Problems beschreibt dabei die Empfindlichkeit der Lösung gegen- überStörungenindenEingangsdaten.DagegenbeschreibtdieStabilitäteinesAlgorithmus,obdieserdie gewünschtenWerteauchunterdemEinflussderunvermeidlichenRundungsfehlergenügendgenauappro- ximiert. Die Einschränkung der beim elastischen Balken auftretenden Minimierungsaufgabe auf einen endlich- dimensionalen Vektorraum von Funktionen führt auf lineare Gleichungssysteme, mit deren numerischer LösungwirunsindenweiterenKapitelndieserVorlesungbeschäftigenwerden.Dabeiunterscheidetman zwischen direkten Verfahren wie dem Gauß-Algorithmus und seinen Varianten sowie iterativen Ver- fahren,diebeiSystemenmitsehrvielenUnbekanntenzumEinsatzkommen.Abschließendwerdennoch iterativeVerfahrenzurnumerischenLösungnichtlinearerGleichungssysteme,wiesiehäufigausmathe- matischenModellenresultieren,behandelt. Kapitel 1 Interpolation von Funktionen Um die für die Auslenkungdes elastischenBalkens in Frage kommenden Funktionen y(x),0 6 x 6 L, mitendlichvielenParameternnäherungsweisedarstellenzukönnen,benötigenwirgeeigneteVektorräume möglichsteinfacherFunktionen.DieallgemeineProblemstellungderInterpolationistdiefolgende: AufeinemreellenIntervall[a,b]seieneineFunktionf: [a,b] RsowieausgewählteStützpunktea 6 → x < x < < x 6 b gegeben.Dazubestimme maneine „einfachere“Funktion p: [a,b] R mit 0 1 n ··· → p(x )=f(x )füri=0,1,...,n. i i 1.1 Polynom-Interpolation Interpolationsaufgabe: Zu paarweise verschiedenenStützpunkten (oder Knoten) x ,x ,...,x [a,b] 0 1 n ∈ bestimmemaneinPolynomp P (vomGradn), n ∈ p(x)=a xn+ +a x+a mit a ,a ,...,a R, n 1 0 0 1 n ··· ∈ welchesfbezüglichderStützpunkteinterpoliert,d.h.p(x )=f :=f(x )füri=0,...,nerfüllt. i i i ZusammengefasstbildendieInterpolationsbedingungenoffenbareinlinearesGleichungssystem(LGS): 1 x ... xn a f 0 0 0 0 1 x ... xn a f ... ...1 ... ...1  ...1 = ...1 ∈Rn+1.      1 x ... xna  f   n n n  n      DurchgeschickteWahleinerBasis{L (x)}n P statt{xi}n erhältmaneinäquivalentesLGSmitder i i=0 ⊂ n i=0 EinheitsmatrixIalsKoeffizientenmatrix.Mansiehtleicht,dassdieLagrange-PolynomeL ,...,L P 0 n n ∈ einesolcheBasisbilden:L istdefiniertalsLösungderspeziellenInterpolationsaufgabe i L (x )=1, L (x )=0 für j=i. i i i j 6 Offenbarerfüllt n x−x L (x):= j i x −x i j j=0,j6=i Y dieseBedingung.DasInterpolationspolynompergibtsichdannunmittelbarzu n p(x)= f L (x). i i i=0 X Damitistdie ExistenzeinerLösungderInterpolationsaufgabe bewiesen.Eindeutigkeit: Seienp,q P n ∈ mitp(x ) = q(x ) = f füri = 0,...,n.DasPolynomp−q P hatdannn+1Nullstellen,alsogilt i i i n ∈ nachdemFundamentalsatzderAlgebrap q. ≡ 4 NumerischeMathematikI Satz1.1. Zun+1paarweiseverschiedenenStützpunktenx mitvorgegebenenStützwertenf ,i=0,...,n,existiert i i genaueinInterpolationspolynomvomGradn. Folgerung1.2. SeiV einFunktionenraum,derP alsUntervektorraumenthält.DannistdiePolynom-InterpolationaufV n bzgl.x ,...,x eineAbbildungP P : V P ,undzwareinelineareProjektion. 0 n ≡ V,x0,...,xn → n Beweis. NachSatz1.1istPAbbildung.Seienf,g V,α Rundp := P(f),q := P(g).Manrechnet ∈ ∈ nach:P(αf)=αp,P(f+g)=p+q(Linearität)undP(P(f))=p=P(f)(Projektionseigenschaft). MitderDarstellungdesInterpolationspolynomsüberdieLagrange-Polynomehabenwiraucheinnumeri- schesVerfahrenzurPolynominterpolation.ZureffizientenAuswertungderFormel n n x−x p(x)= f j i x −x i j i=0 j=0,j6=i X Y fürbeliebigex [a,b]sindjedocheinigeÜberlegungenangebracht.WirleitenimfolgendeneinenAlgo- ∈ rithmuszureffizientenAuswertungdieserFormelher. DazuführenwirdasKnotenpolynomω (x):=(x−x ) (x−x )einundsehen,dassdieLagrange- n+1 0 n ··· PolynomedieIdentität n ω (x) 1 L (x)= n+1 i x−x x −x i i j j=0,j6=i Y erfüllen.MitdenGewichten n 1 w := , i=0,...,n, i x −x i j j=0,j6=i Y lautetdieresultierendeFormelfürdieAuswertungdesInterpolationspolynomsaneinerStellex: n n ω (x) f w p(x)= f n+1 w =ω (x) i i . i x−x i n+1 x−x i i i=0 i=0 X X Algorithmus(BerechnungderGewichtew fürdiePolynominterpolation). i fork=0tondo fori=0tok−1do w =w (x −x ) i i i k endfor w =1 k fori=0tok−1do w =w (x −x ) k k k i endfor endfor fork=0tondo w =1/w k k endfor DerAufwandfür jede AuswertunganeinemPunkt xbeiebenfalls variablenStützwertenf beträgt nach i anfänglicherBerechnungdern+1Gewichtew :2n+1Additionenbzw.Subtraktionen,3n+3Multiplika- i tionenbzw.Divisionen.FürgegebeneStützwertef reduziertsichderAufwandumn+1Multiplikationen, i wennmandieProduktef w vorabberechnet. i i Satz1.3(FehlerdarstellungzurPolynom-Interpolation). DieFunktionfsei(n+1)-malstetigdifferenzierbar,f Cn+1[a,b].Danngibteszujedemx [a,b]ein ∈ ∈ ξ (a,b)so,dassfürdasInterpolationspolynomp(x)bezüglichderKnotenx ,x ,...,x gilt: 0 1 n ∈ f(n+1)(ξ) f(x)−p(x)= ω (x). (n+1)! n+1 1.InterpolationvonFunktionen 5 1 0.5 0 -0.5 -1 -1 -0.5 0 0.5 1 Abbildung1.1:Tschebyschow-PolynomeT ,...,T 0 5 Beweis. Wirhaltenx¯ [a,b]festundbetrachtendenFallx¯ =x ,i=0,...,n.Setze i ∈ 6 f(x¯)−p(x¯) F(x):=f(x)−p(x)− ω (x). ω (x¯) n+1 n+1 Fbesitztmindestensn+2Nullstellenin[a,b]: f(x¯)−p(x¯) F(x )=f(x )−p(x )− ω (x )=0, i=0,...,n, i i i ω (x¯) n+1 i n+1 =0 f(x¯)−p(x¯) F(x¯)=f(x¯)−p(x¯)− ω| {(xz¯)=}0. ω (x¯) n+1 n+1 NachdemSatzvonRollebesitztF′(mindestens)n+1Nullstellenin(a,b),...,F(n+1)(mindestens)eine Nullstelleξin(a,b).Esfolgt f(x¯)−p(x¯) 0=F(n+1)(ξ)=f(n+1)(ξ)−p(n+1)(ξ)− ω(n+1)(ξ) ω (x¯) n+1 n+1 =0 =(n+1)! = f(x¯)−p(x¯)|= {f(zn+1})(ξ)ω (x¯) | {z } ⇒ (n+1)! n+1 Bemerkung. AusSatz1.3ergibtsichunmittelbar: 1 f(n+1) |f(x)−p(x)|6 (n+1)!ξ∈m[aax,b]|f(n+1)(ξ)|·|ωn+1(x)|= k(n+1k)!∞|ωn+1(x)|. Definition(Tschebyschow-Polynome1.Art). DieTschebyschow-Polynome1.Art(Abb.1.1)sindrekursivdefiniertdurch: T (t):=1, T (t):=t, T (t):=2tT (t)−T (t), k=1,2,... 0 1 k+1 k k−1 Lemma1.4(EigenschaftenderTschebyschow-Polynome1.Art). (a) Fürt [−1,1]giltT (t)=cos(karccost). k ∈ (b) Fürk>1giltT (t)=2k−1tk+c tk−1+ +c t+c mitc ,...,c R,d.h.,derführende k k−1 1 0 k−1 0 ··· ∈ KoeffizientvonT (t)ist2k−1. k 6 NumerischeMathematikI (c) Esgilt|T (t)|61fürt [−1,1]. k ∈ Beweis. (a)Fürk=0undk=1istdieFormeloffenbarrichtig.Weitergilt cos((k 1)arccost)=cos(karccost)cos(arccost) sin(karccost)sin(arccost) ± ∓ = cos((k+1)arccost)+cos((k−1)arccost)=2tcos(karccost), ⇒ d.h.,dieRekursionsvorschriftfürdieTschebyschow-Polynome.(b)folgtdirektausderRekursionsformel, (c)folgtaus(a). Beispiel. NullstellenvonT (x)alsStützpunktefürdiePolynom-Interpolationauf[−1,1]: n+1 1 2j+1 T (x)=0 (n+1)arccosx= j+ π, j Z x=cos π , j Z. n+1 ⇐⇒ 2 ∈ ⇐⇒ 2(n+1) ∈ (cid:18) (cid:19) (cid:18) (cid:19) Umsortierung:j=n−i 2(n−i)+1 2i+1 x =cos π =−cos π , i=0,1,...,n. i 2(n+1) 2(n+1) (cid:18) (cid:19) (cid:18) (cid:19) DaszugehörigeKnotenpolynomist(vgl.Lemma1.4(b)) 1 1 ω (x)= T (x) undnachLemma1.4(c)gilt max |ω (x)|= . n+1 2n n+1 x∈[−1,1] n+1 2n ZumVergleichbeiderVerwendunggleichverteilterStützpunkte n i i x¯ =−1+2 , ω¯ (x)= x+1−2 . i n n+1 n i=0(cid:18) (cid:19) Y n 5 10 20 max |ω¯ (x)| 6.9 10−2 8.5 10−3 2.3 10−4 x∈[−1,1] n+1 · · · max |ω (x)| 3.1 10−2 9.8 10−4 9.5 10−7 x∈[−1,1] n+1 · · · Beispiel. WirbetrachtendieInterpolationderFunktionf C∞[−1,1]mit ∈ 0 fürx [−1,0], f(x):= ∈ (e−1/x fürx (0,1]. ∈ DieAbbildungen1.2und1.3zeigendeutlich,dassdieDifferenzzwischendemInterpolationspolynomund derzuinterpolierendenFunktionfbeiVerwendungderTschebyschow-Knotenerheblichkleinerist. DassdieNullstellenderTschebyschow-PolynomebezüglichderAbschätzungausSatz1.3dieoptimalen KnotenfürdiePolynominterpolationsind,istdieAussagedesfolgendenSatzes. Satz1.5. FürjedesPolynomq derForm n+1 q (x)=(x−x )(x−x ) (x−x ) n+1 0 1 n ··· gilt: 1 1 max |q (x)|> = max |T (x)| . x∈[−1,1] n+1 2n 2n x∈[−1,1] n+1 (cid:18) (cid:19) 1.InterpolationvonFunktionen 7 0.35 f(x) gleichm. 0.3 Tscheb. 0.25 0.2 0.15 0.1 0.05 0 -1 -0.5 0 0.5 1 Abbildung1.2:PolynominterpolationvomGrad5 0.4 f(x) gleichm. 0.35 Tscheb. 0.3 0.25 0.2 0.15 0.1 0.05 0 -1 -0.5 0 0.5 1 Abbildung1.3:PolynominterpolationvomGrad10 Beweis(Widerspruch). Angenommen,esgibteinPolynomq derForm n+1 q (x)=(x−x )(x−x ) (x−x ), n+1 0 1 n ··· sodass|q (x)| < 2−n fürallex [−1,1]gilt.Dannistp(x) := q (x)−2−nT (x) P ,also n+1 n+1 n+1 n ∈ ∈ einPolynomvomGradn,fürdasandenPunktenξ =cos(jπ/(n+1)),j=0,...,n+1gilt: j p(ξ )=q (ξ )−2−ncos((n+1)arccosξ ) j n+1 j j =q (ξ )−2−ncos(jπ)=q (ξ )−2−n(−1)j, j=0,...,n+1. n+1 j n+1 j 8 NumerischeMathematikI Wirhabensomit <+2−n−2−n =0 fallsjgerade, p(ξ )=q (ξ )−2−nT (ξ ) j n+1 j n+1 j (>−2−n+2−n =0 fallsjungerade. Diesbedeutet,dassp injedemderIntervalle (ξ ,ξ ),j = 0,...,neineNullstelle besitzt.Somit muss j j+1 p 0geltenimWiderspruchzurAnnahme. ≡ 1.2 Spline-Interpolation BeiderPolynominterpolationwurdedieApproximationsgütedurchErhöhungdesPolynomgradesverbes- sert.DieseVorgehensweisehatjedochfolgendeNachteile: dieAbschätzungfürdenInterpolationsfehlerhängtvonAbleitungenhöhererOrdnungab, • dieInterpolationspolynomeneigenzustarkenOszillationen. • Grundidee der Spline-Interpolation: berechne zwischen den Knoten separate Polynome mit identischem (und relativ kleinem) Polynomgrad, so dass in den Knoten die Funktion interpoliert wird und zusätzlich geeigneteÜbergangsbedingungenerfülltsind. Definition(Spline). Zun+1Stützpunktena = x < < x = bistdurchT := {[x ,x ],...,[x ,x ]}eineZerlegung 0 n 0 1 n−1 n ··· desIntervalls[a,b]gegeben.EinSplinevomGradkbezüglichTisteineFunktions Ck−1[a,b],dieauf ∈ jedemIntervall[x ,x ],i = 1,...,n,miteinemPolynoms P übereinstimmt.Wirbezeichnenden i−1 i i k ∈ RaumderSplinesvomGradkbezüglichTmitS . k,T Beispiel(k=1,linearerSpline,dimS =n+1). 1,T (cid:8)t (cid:8)s(cid:8)1 (cid:8)`t ``s2``t s3 tXsX4 `t ``s5``t a=x x x x x x =b 0 1 2 3 4 5 Lemma1.6. FürdieDimensiondesSpline-RaumesS gilt:dimS 6n+k. k,T k,T Beweis. In[x ,x ]betrachtenwireinbeliebigesPolynoms P ,alsok+1freieParameter.Fürs (x), 0 1 1 k i i = 2,...,n sind jeweils die k Anfangswerte s (x ),s′(x∈ ),...,s(k−1)(x ) durch die Bedingung i i−1 i i−1 i i−1 s C(k−1)[a,b] festgelegt. Somit bleibt für jedes s , i = 2,...,n, ein freier Parameter, also insgesamt i ∈ (k+1)+(n−1)=n+k. WirwendenunsdemSpezialfallkubischerSplines,d.h.k=3zu.NachLemma1.6giltdimS 6n+3. 3,T Wirgebenn+3Bedingungenvorundzeigenspäter,dassdieseAufgabeeineLösungbesitzt.Darausfolgt danndimS =n+3. 3,T Nebendenn+1Interpolationsbedingungens(x )=f ,i=0,...,n,sindfolgendeMöglichkeitenfürdie i i restlichenzweiBedingungenüblich: (A) s′(a)=f′,s′(b)=f′ zuvorgegebenenWertenf′,f′ R(vollständigerSpline), 0 n 0 n ∈ (B) s′′(a)=0,s′′(b)=0(natürlicherSpline), (C) s′(a)=s′(b),s′′(a)=s′′(b)unterderZusatzvoraussetzungf =f (periodischerSpline). 0 n