Persönliche Website – Services & Ressourcen | ETH …ricklis/blog/download... · Web...

a b

1 2 3 4

a b

1 2 3 4

Zusammenfassung Numerische Methoden Stefan Rickli http://blogs.ethz.ch/ricklis

1 Numerische Quadratur

Ziel: Berechne Wert von ∫a

b

f ( x )⋅dx , wenn es nicht analytisch lösbar ist

⇒∫a

b

f (x ) ⋅dx=∫a

b

∑i=0

n−1

αi bi ( x )=∑i=1

n

ωi f (c i)

ωi: Gewichte

c i: Knoten

1.1 Transformation auf eine ReferenzformelWollen das Integral mit allgemeinen Grenzen auf ein Referenzintervall bringen, um standardisierte Formeln verwenden zu können.

[a ,b ]→ [−1,1 ]

∫a

b

f ( t ) ⋅dt= 12

(b−a )∫−1

1

f (τ ) ⋅dτ ≈ 12

(b−a )∑j=1

n

ω j f ( c j )=∑j=1

n

ω j f (c j )

c j=b−a

2c j+

a+b2

,ω j=12

(b−a ) ωj

1.2 Newton-CotesKnoten sind äquidistant auf [a ,b ] verteilt

1.2.1 Mittelpunktregel

QM ( f ;a , b )=(b−a ) f ( a+b2 )Fehler:

E=|∫ab

f ( x ) ⋅dx−QM ( f ; a , b )|= 124

(b−a )3|f (2 ) (ξ )|Konvergenzordnung: 2

1.2.2 Trapez-Regel

QT ( f ; a , b )=b−a2 ( f (a )+ f (b ) )

Fehler:

E= 112

(b−a )3|f ( 2) (ξ )|

Konvergenzordnung: 2

1.2.3 Simpson Regel

QS= b−a6 (f (a )+4 f ( a+b2 )+ f (b ))

Fehler:

E= 190 ( b−a

2 )5

|f ( 4 ) (ξ )|Konvergenzordnung: 4

1.3 Summierte QuadraturProblem: [a ,b ] ist gross → Näherung schlechtIdee:

Erhöhe Genauigkeit durch Polynome höheren Grades → negative Gewichte ab Ordnung 8 ⇒ Fehler explodiert

Verwende summierte QuadraturUnterteile [a ,b ] in N Teilintervalle, wende auf jedes die

Quadratur an

∫a

b

f ( x )⋅dx=∑i=0

N−1

∫xi

x i+1

f ( x )⋅dx≈∑i=0

N−1

Q (f ; x i , x i+1 )

h=b−aN

, xa=a , x i=x0+i⋅h , xN=b

1.3.1 summierte Mittelpunkt-Regel

∫a

b

f ( x )⋅dx ≈∑i=0

N−1

hf ( xi+x i+1

2 )Fehler:

E (h )≤ h2

24(b−a ) max

x∈ [ a ,b ]{|f ' ' ( x )|}

1.3.2 summierte Trapez-Regel

∫a

b

f ( x )⋅dx ≈ h2 ∑

i=0

N−1

( f (x i )+ f (xi+1 ))=h2 [ f (a )+2∑

i=1

N−1

f (xi )+ f (b )]Fehler:

E (h )≤ h2

12(b−a ) max

x∈ [a,b ]{|f ' ' (x )|}

1.3.3 summierte Simpson-Regel

∫a

b

f ( x )⋅dx∧≈ h6 ∑

i=0

N−1 ( f (x i )+4 f ( x i+x i+1

2 )+ f (xi+1 ))¿= h

6 [ f (a )+2∑i=1

N−1

f (xi )+4∑i=1

N

f ( x i−1+x i

2 )+ f (b )]Fehler:

E (h )≤ h4

2880(b−a ) max

x∈ [a , b]{|f (4 ) (x )|}

zuletzt geändert: 03.01.2017 19:42, Version 20 1 / 16


1.3.4 Adaptive QuadraturUm eine gewisse Genauigkeit zu erreichen, muss die Schrittweite klein genug gewählt werden. Dies zieht einen hohen Rechenaufwand nach sich, vor allem wenn die zu integrierende Funktion aufwändig auszuwerten ist. Deshalb versucht man bei vorgegebener Genauigkeit die Schrittweite möglichst gross zu wählen.

1) Werte Quadratur für das nächste Teilintervall mit Trapez- (Ordnung 2) und Simpson-Regel (Ordnung 4) aus.

2) Bilde Differenz zwischen den beiden Resultaten.3) Falls die Differenz einen bestimmten Toleranzbetrag überschreitet,

halbiere das Intervall4) Führe Schritte 1-3) noch einmal auf den neuen Teilintervallen aus.5) Falls die Toleranz eingehalten wird, kann man den Wert der

Simpson-Regel für das Teilintervall verwenden.

Beispiel für adaptive Quadratur mit Hilfe der Trapez- und Simpsonregel:f = @(x) exp(6*sin(2*pi*x)); %Integrandxstart = 0;xend = 0.9; %Intervallrefval = integral(f,xstart,xend); %Referenzwerttol = 1e-3; %Toleranz fuer Adaptivitaethtoosmall = false; %gehe aus while-loop fuer h zu

kleinfstart = f(xstart); FEND = f(xend); %Start- und Endwert von fapprintegral = 0; %Approximation des Integraleswhile xstart < xend

h = xend - xstart; %Versuche mit einem grossen Schritt

fend = FEND; %Werf von f am Endpunkt und in der Mitte des Intervalls

fmiddle = f(xstart + h/2);%Berechne approximierte Loesung mit Trapez- und Simpsonregelhfse = h*(fstart + fend);trapez = hfse/2;simpson = (hfse+4*h*fmiddle)/6;err = abs(simpson - trapez); %absoluten Fehler%wenn Fehler zu grosswhile err > tol

h = h/2;%halbiere Schrittweite%when h zu klein ist, dann unterbreche die Rechnung. Wahrscheinlich ist etwas falsch implementiert worden.if h < 1e-8

htoosmall = true; break

endfend = fmiddle; %Werf von f am Endpunkt und in

der Mitte des Intervallsfmiddle = f(xstart + h/2);hfse = h*(fstart + fend); %Berechne

approximierte Loesung mit Trapez- und Simpsonregel

trapez = hfse/2;simpson = (hfse+4*h*fmiddle)/6;err = abs(simpson - trapez); %absoluter Fehler

endapprintegral = apprintegral + simpson; %Addiere

die Approximationxstart = xstart + h; %Update xstart und fstartfstart = fend;%when h zu klein ist, dann unterbreche die Rechnung, Wahrscheinlich ist etwas falsch implementiert worden.if htoosmall

apprintegral = NaN; fprintf('h ist zu klein \n')break

endendfprintf('Approximation = %2.4f, Referenzwert = %2.4f\n',apprintegral,refval)

1.4 Gauss-QuadraturGauss-Quadraturformeln haben nicht äquidistante Stützstellen. Mit n Knoten wird Ordnung 2n erreicht. D.h. ein Polynom vom Grad

2n−1 wird exakt integriert.

Beispiel:Wollen 2 Knoten und Ordnung 4, d.h. 1 , x , x2 , x3 exakt integrieren

¿∫−1

1

1⋅dx=2=1w1+1w2

¿∫−1

1

x ⋅ dx=0=c1w1+c2w2

¿∫−1

1

x2⋅ dx=23=c1

2w1+c22w2

¿∫−1

1

x3 ⋅dx=0=c13w1+c2

3w2

Löse die 4 Gl. nach w1 ,…,c2 auf

⇒w1=w2=1 und c1=√ 13, c2=−√ 1

3

∫−1

1

f ( x ) ⋅dx≈1 ⋅ f (√ 13 )+1⋅ f (−√ 1

3 )Anzahl Knoten n Knoten c i Gewichte ωi

1 0 2

2 ±√ 13

1

3

0

±√ 35

8959

(Gauss-Legendre-Quadraturformeln haben immer positive Gewichte)

1.5 OrdnungEine Quadraturformel hat die Ordnung n+1, wenn damit Polynome der

Ordnung n exakt integriert werden.

Theorem: Die Ordnung einer symmetrischen Quadraturformel ist gerade.

Ordnung einer Quadratur mathematisch bestimmen: Monome mit steigender Ordnung in Quadraturformel einsetzen und mit analytischer Lösung des entsprechenden Integrals vergleichen, bis zum ersten Mismatch. Aufgrund der Linearität der Quadratur kann argumentiert werden, dass auch Polynome bis zu diesem Grad exakt integriert werden.



1.6 KonvergenzFür eine zusammengesetzte Quadraturformel Q, die auf einer lokalen Quadraturformel der Ordnung p∈N basiert, gibt es ein Konstante C>0, sodass

E (h )=|∫ab

f ( t )⋅ dt−Q ( f ;a , b )|≤C hp‖f ( p )‖L∞ ( [a ,b ] )∀ f ∈C p ( [a ,b ] )

Genauigkeit vergrössern: Ordnung der Formel erhöhen: z.B. MP → Simpson kleinere Teilintervalle

Für volle Ordnung muss f hinreichend glatt sein!

1.6.1 algebraische KonvergenzBei der algebraischen Konvergenz verhält sich der Fehler wie

E (h )=O (hp )Grafische Bestimmung der Konvergenzordnung:benütze log-log Plot, vergleiche mit h ,h2 , h3 ,… , verwende least-squares Fit:fit = polyfit(log(h),log(E(h)),1);fprintf('Konvergenzrate = %1.2f\n',fit(1));

Betrachte die Fehler von E (h ) und E( h2 ):

E (h )E (h/2 )

= Chp

C (h/2 )p=2p⇔ p=

ln (E (h ) /E (h /2 ) )ln (2 )

Beispiel: Summierte Quadratur

1.6.2 exponentielle KonvergenzBei der exponentiellen Konvergenz verhält sich der Fehler wie

E (n )=O (qn ) ,0<q<1

Grafische Bestimmung der Konvergenzordnung:benütze semi-log Plot, benütze least-squares Fit:fit = polyfit(n,log(E(n)),1)fprintf('Konvergenzrate = %1.2f\n',exp(fit(1)));

Betrachte die Fehler von E (m ) und E (n ):E (m )E (n )

=Cqm

Cqn =qm−n⇔ ln (q )=ln (E (m )/E (n ) )

m−n

Beispiel: Gauss-Legendre Quadratur

2 gewöhnliche Differentialgleichungen2.1 ODE und AWPODE 1 Ordnung:

y ( t )=f (t , y ( t ) )ODE n-ter Ordnung:

y (n) (t )=f (t , y (t ) , y ' ( t ) ,…, y ( n−1) (t ) )Umwandlung in ein ODE-System 1. Ordnung

z=(z1

z2

⋮zn

)=(yy⋮

y (n−1))⇒ z ( t )=(yy ( 2)

⋮y (n))=¿

Anfangswertproblem (AWP):

ODE n-ter Ordnung und y (t 0 )= y0 , y' (t 0 )= y 0

' ,… gegeben.

Autonomitätf hängt nicht explizit von Zeit ab, d.h.

f ( t , y (t ) )=f ( y (t ) )

2.2 Existenz und Eindeutigkeit der LösungLipschitz-Stetigkeit:

Eine Funktion f :R→R heisst lipschitz-stetig, falls

∃L>0 :|f (x2 )−f (x1 )|≤ L|x2−x1|

2.2.1 Satz von Piccard-Lindelöf:Falls

f stetig in (t , y ) (in t‼‼)

f lipschitz-stetig in y auf [ t 0 ,t 0+δ ]×D in Umgebung des

Anfangswerts y0⇒ Es existiert eine eindeutige Lösung y ( t ) von

y ( t )=f (t , y ( t ) ) , y (t 0 )= y0 für zumindest eine kurze Zeitspanne

[ t 0 ,t 0+ε ] , ε>0

Bemerkungen:Lipschitz-Stetigkeit von f in y ist WICHTIG

f stetig, f ' beschränkt ⇒ lipschitz

f stetig diff’bar ⇒ lipschitz

f nicht stetig ⇒ f nicht lipschitz Kegel malen (wo man sieht, ob Steigungstangente geschnitten wird

(⇒ nicht lipschitz) oder nicht (⇒ lipschitz))

Für globale Lösungen muss noch angenommen werden, dass f beschränkt bleibt.



3 LöserGegeben: AWP y ( t )=f (t , y ( t ) ) , y (t 0 )= y0

3.1 Runge-Kutta MethodenAnalytische Lösung des AWPs:

y (t1 )= y0+∫t 0

t 1

~f ( τ ,~y (τ ) ) ⋅dτ

Berechne das Integral numerisch mit Quadratur:

∫0

1

g ( x )⋅ dx=∑i=1

s

ωi g (c i )(Zuerst muss noch auf ein Referenzintervall transformiert werden).

∫t0

t1

f (τ , y ( τ ) ) ⋅dτ =h≔t 1−t0 ¿∫

0

h

f (τ+t 0 , y (τ+t 0 )) ⋅dτ

≈∧∑ihωi f (t 0+c ih , y (t 0+ c ih ))

⇒ y (t 1 )≈ y0+h∑iωi f (t 0+ci h , y (t 0+ cih ) )

Definition Runge-Kutta Verfahren:

¿bi , aij∈R ,c j=∑j=1

s

aij

¿k i=f (t n+ci h , yn+h∑j=1

s

aijk j)¿ yn+1= y n+h∑

i=1

s

bi k i}s-stufiges R-K Verfahren

Butcher-Tableau:

c A c1

c2

⋮cs

a11 … a1 s

a21 ¿ ¿⋱¿ ⋮ ¿as1¿…¿ass¿

bT b1 … bs

A∈R s×s: Runge-Kutta-Matrix

b∈R s: Gewichte

c∈R s: Knoten

Es gilt ∑i=1

s

bi=1 und

c i=∑j=1

s

aij

3.2 Konvergenz von RK-Verfahren / Diskretisierungsfehler3.2.1 globaler DiskretisierungsfehlerSei y : [ t0 , T ]→Rn Lösung eines AWP. Dann ist

eh ( tn )= y (tn )− ynapprox

der Fehler zum Zeitpunkt t n, der sich in n Schritten akkumuliert hat.

Ein R-K-Verfahren heisst konvergent der Ordnung p, falls

maxk=0 ,…, n

‖eh (t k )‖=O (hp ) , h→0



3.2.2 lokaler DiskretisierungsfehlerDer lokale Diskretisierungsfehler ist der Fehler, den man in einem Schritt macht, weshalb er auch Ein-Schritt-Fehler genannt wird.

~yn+1∧≔ y (t n )+h∑i=1

s

bi k i

ρn+1∧¿ y (tn+1 )−~yn+1

~yn+1 basiert auf exaktem Wert zum Zeitpunkt t n

Der Abschneidefehler oder Konsistenzfehler ist gegeben durch

τ n+1=ρn+ 1

h=

y (tn+1 )− y (t n )h

−h∑i=1

s

bi k i→ benütze Taylorentwicklung

Ein R-K-Verfahren hat die Konsistenzordnung p, falls

maxk=0 , …, n

{‖τ k‖}=O (hp ) , h≤h0

Es gilt, falls f lipschitz:Konsistenz der Ordnung p ⇒ Konvergenz der Ordnung p

Konsistenz + Stabilität ⇒ Konvergenz

Verfahren Konsistenz-OrdnungExpl. Euler 1Expl. Trapez 2Expl. Mittelpunkt 2Klassisches RK4 4Kutta’s 3/8 4

3.3 explizite Löser3.3.1 Expliziter Euler (Ordn. 1)Verwende Taylorapproximation 1. Ordnung: y (t+h )≈ y (t )+h y (t)Verwende Definition der DGL: y ( t )+h y ( t )= y (t )+h f (t , y (t ) )Iterationsvorschrift:

yn+1= yn+h f (t n , yn )

3.3.2 Explizite Trapezregel (Ordn. 2) Verwende expl. Euer in

yn+1= yn+h( 12f ( tn , yn )+ 1

2f (t n+1 , yn+1 ))

für Berechnung vonyn+1 der rechten Seite.

k1∧¿ f (t n , yn )k 2∧¿ f (t n+h , yn+hk1 )

yn+1∧¿ yn+h2 (k 1+k2 )

3.3.3 Explizite Mittelpunktregel (Ordn. 2)

k1∧¿ f (t n+c1h , yn+ha11k1 )∧¿=f (t n , yn )k 2∧¿ f (t n+c2h , yn+ha21 k1+ha22 k2 )∧¿=f ( tn+h /2 , yn+h/2k1 )

yn+1∧¿ yn+h⋅ k2

3.3.4 weitere Verfahren

klassisches RK4 Kutta’s 3/8-Regel0 0 0 0 012

12

0 0 0

12

012

0 0

1 0 0 1 0

0 0 0 0 013

13

0 0 0

23

−13

1 0 0

1 1 −1 1 0


0 01

0 0 01 1 0

12

12

0 0 012

12

0

0 1


16

26

26

16

18

38

38

18

3.4 implizite LöserRechte Seite der Iterationsvorschrift hängt ebenfalls von yk +1 ab.⇒ führt auf ein implizites Gleichungssystem, i.d.R nichtlinear

yn+ 1= yn+h f (t n+1 , yn+1 )z.B. wenn y=sin ( y )⇒ yn+1= yn+h sin ( yn+1)Lösen durch Newtonverfahren

Suche Nullstelle von φ ( y )= y0+h sin ( y1 )− y1

Verfahren Schritt Ordn.

Impl. Euler yn+ 1= yn+h f (t n+1 , yn+1 ) 1

Impl.Trapezregel yn+ 1= yn+

h2 [ f ( tn , yn )+f (t n+1 , yn+1 ) ] 2

Impl.Mittelpunktregel yn+ 1= yn+h f ( 1

2 ( tn+t n+1 ) , 12 ( yn+ yn+1 )) 2

3.5 Semi-Implizite LöserBasieren auf der Grundlage von impliziten Verfahren. Anstatt aber das Gleichungssystem zu lösen, wird nur ein Newton-Schritt durchgeführt.

Implizites Verfahren in ein Semi-Implizites Verfahren überführen:1. yn+1=… umformen in ein Nullstellenproblem

ϕ ( yn , yn+1 )=02. Newton-Schritt formulieren:

yn+1(1 ) = yn+1

(0 ) −ϕ ( yn , yn+1

( 0) )∂

∂ yn+1ϕ ( y n , yn+1

(0 ) )3. Benütze als Startwert yn für yn+ 1

( 0)

⇒ yn+1( 1) = yn−

ϕ ( yn , yn )ϕ ' ( yn , yn )



4 Newton-Verfahren4.1 FixpunktiterationSchreibe F ( x )=0 um in Fixpunktproblem ϕ ( x )=xDies liefert die Iterationsvorschrift x (k +1)=ϕ (x ( k ) )ϕ ( x ) ist in der Regel nicht eindeutig und die Konvergenz unterscheidet sich je nach gewählter Iterationsvorschrift.

function x = fixIt(F, x0, tol, kmax)x(1) = x0;for i=1:kmax

x(i+1) = F(x(i));if abs(x(i+1) - x(i)) < tol

break; end; end;

Konvergiert eine Iteration lokal gegen den Fixpunkt x¿?

|ϕ' ( x¿ )|<1∧⇒ lokale Konvergenz|ϕ' ( x¿ )|≥1∧⇒konvergiert nicht

Def: Eine Iterationsvorschrift heisst konsistent, falls

f ( x¿ )=0⇔ϕ (x¿)=x¿

gilt.

Def: ϕ ( x ) heisst Kontraktion auf D⊂Rn, falls ∃L mit 0<L<1, sodass

‖ϕ ( x )−ϕ ( y )‖≤ L‖x− y‖∀ x , y∈D

4.1.1 Banachscher FixpunktsatzFalls ϕ eine Kontraktion ist und eine abgeschlossene Menge M⊂D auf sich selbst abbildet, dann gilt

i) ϕ hat genau einen Fixpunkt in Mii) für jeden Startwert x (0 )∈M gilt: x (k )=ϕ (x (k−1) )

konvergiert gegen x¿

iii) Fehlerschätzung: ‖x ( k )−x¿‖≤ Lk

1−L‖x ( 1)−x (0 )‖

4.2 Startwertx (0 ) muss nahe genug an x¿ liegen

Falls x (0 ) zu weit weg liegt oder x (k +1) weiter weg von x¿ liegt als x (k ):benütze das gedämpfte Newtonverfahren.Dämpfe Schnittpunkt der Tangente mit x-Achse, falls dieser weiter entfernt ist als vorher.



4.3 Newton-Verfahren 1DIdee: Starte mit x (0 ) und berechne x (k +1) aus x (k )

- berechne Nullstelle der Tangente- Nullstelle ist neues x (k )

φ ( x )≔F ( x )Taylorentwicklung:

F ( x )≈F (x (k ))+F ' ( x( k ) ) (x−x (k ) )Nullstelle von F ( x ) :

F (x ( k+1 ))=F (x (k ) )+F' (x (k )) (x ( k+1)−x( k ) )=! 0

x (k +1)=x ( k )−F (x ( k ) )F' (x (k ))

Voraussetzungen für erfolgreiche Iteration:

- F ' (x (k ) )≠0- x (0 ) muss nahe genug an der gesuchten Nullstelle liegen

Algorithmus:function x=newton(x0,F,DF,nMax,TOL)x = x0;Fx_new = F(x);for ii = 1:nMax %Newton iterationen

Fx = Fx_new; %Wert von F undDFx = DF(x); %DF in alter Iterierteif abs(DFx)<1e2*eps %Bedingung: DF nicht Null

fprintf('DFx ist zu klein\n')assert(false);

enddelta_x = - DFx\Fx;%Newton-Updatex = x + delta_x; %Neue IterierteFx_new = F(x); %Wert von F in neuer Iterierteif abs(DFx\Fx_new) < TOL; %wenn update zu klein

break % stop and exitend

endend

4.4 gedämpftes Newton-Verfahren

x (k +1)=x ( k )−λ (k ) F (x ( k ) )F ' ( x( k ) )

, λ(k )>0



4.5 Newton-Verfahren NDSei F :Rn→Rn zweimal stetig diff’bar

Wir suchen x¿ mit F (x¿)=0

Wie im 1D-Fall: Ersetze F durch ein linearisiertes ~F~F( k ) ( x )=F (x (k ) )+DF (x (k )) (x−x ( k ) )

bestimme x (k +1) :F ( x( k ) )+DF (x ( k ) ) (x (k+1 )−x (k ) )=0⇔DF ( x( k ) ) (x (k +1)−x ( k ) )⏟

Δx( k )

=−F (x (k ) )

1) Stelle Jakobideterminante auf

2) löse DF (x (k ) ) Δ x( k )=−F (x (k ) )3) Update x (k +1)=x ( k )+Δ x (k )

in einem Schritt:

x (k +1)=x ( k )−DF (x (k ) )−1F (x (k ))ABER: es ist teuer, Inverse zu berechnen, weshalb in Praxis voriges Vorgehen mit Lösen des Gleichungssystems zum Einsatz kommt.

4.6 Abbruchkriterien„Wann sind wir nahe genug an x¿?“

Möglichkeit 0:Gar nicht genug… Abbruch aufgrund Überschreitung einer gewissen Anzahl Iterationen (wahrscheinlich Divergenz).

Möglichkeit 1:

‖F ( x( k+1 ) )‖≤TOLProblem: Bei gegebener Toleranz würde man bei einem Problem mit z.B.~F=2F tendenziell mehr Iterationsschritte machen, obwohl die Nullstelle für beide Probleme gleich ist, weil im k-ten Iterationsschritt der „Fehler“ zur Toleranz noch doppelt so gross ist.Man sagt, dass dieser Ansatz nicht affin-invariant ist.

Möglichkeit 2:Skaliere mit der Jacobi-Matrix, um die Iterierte zu normieren:

‖DF (x ( k+1) )−1F (x ( k+1) )‖≤TOL

Dies ist sehr aufwendig, aber wir kennen schon DF (x (k ) )−1

⇒‖DF ( x( k ) )−1F (x ( k+1) )‖≤TOL

4.7 Konvergenz OB ein Verfahren konvergiert, hängt nicht von der Wahl der Norm

ab. Die Wahl der Norm hat Einfluss darauf, in welche

Konvergenzordnung ein Verfahren eingestuft wird.

1-Norm ‖v‖1=|v1|+|v2|+|v3|2-Norm (Euklidische Norm) ‖v‖2=√v1

2+v22+v3

2

∞-Norm / Max-Norm ‖v‖∞=max {|v1|,|v2|,|v3|}4.7.1 lineare KonvergenzEine Folge x (k ) heisst linear konvergent bezüglich x¿, falls

∃0<L<1:‖x( k+1 )−x¿‖≤ L‖x( k )−x¿‖∀ k≥ k 0

4.7.2 quadratische KonvergenzEine Folge x (k ) heisst quadratisch konvergent bezüglich x¿, falls

∃0<C :‖x ( k+1)−x¿‖≤C‖x (k )−x¿‖2∀ k ≥k0

4.7.3 Konvergenz des 1D-Newton-VerfahrenWie schnell konvergiert x (k ) gegen x¿?

Annahmen:i) wähle x (0 ) nahe an x¿

ii) F ' ( x ) habe eine beschränkte Inverse in einer Umgebung

von x¿

iii) sei F zweimal stetig differenzierbar und |F ' ' ( x )| sei

beschränkt in einer Umgebung von x¿

ii) und iii) sind ausreichend für die Aussage, dass F eine Kontraktion ist ⇒ können Banach’schen Fixpunktsatz anwenden

⇒ Die Newtoniterierten x (k ) konvergieren quadratisch gegen x¿.

numerisch lässt sich die Konvergenz wie folgt berechnen:

p≈ ln (e ( k+1)/e (k ) )ln (e( k ) /e (k−1 ))

, e (k )=|x (k )−x¿|Man kann beobachten, dass sich die Anzahl der exakten Dezimalstellen mit jedem Schritt (bis zur Maschinengenauigkeit) verdoppelt.

4.7.4 Konvergenz des ND-Newton-Verfahren- quadratisch

- wenn bei Newtonkorrektur DF (x (k−1 )) verwendet wird:

DF (x (k−1 )) Δ x ( k )=−F (x (k ))⇒ Konvergenz wird linear

(LR-Zerlegung von DF (x (k−1 )) kann wiederverwertet werden)

4.8 Quasi-Newton-VerfahrenKommen bevorzugt in Optimierungsproblemen zum Einsatz, wenn in

∇ f ( x¿)=0x¿ gesucht wird. Die Jacobi-Matrix des Newton-Verfahrens entspricht dann der Hesse-Matrix der Funktion f. Das Aufstellen dieser ist sehr aufwendig, weshalb geeignete Approximationen in diesen Quasi-Newton genannten Verfahren verwendet werden.



5 SchrittweitensteuerungGegeben sei eine Näherung yn zum Zeitpunkt t n sowie ein Vorschlag

für die Schrittweite hn+1=t n+1 – t n für die nächste Näherung yn+1 :

1. Berechne eine lokale Toleranz ¿ ln+1 für den Ein-Schritt-Fehler auf

dem Zeitintervall [ t n ,t n+1 ]2. Mache einen Zeitschritt mit dieser Schrittweite; schätze den Ein-

Schritt-Fehler:a. Ist der Fehler kleiner als ¿ ln+1 , akzeptiere den Schrittb. Ist der Fehler größer, wiederhole den Schritt mit einer

kleineren Schrittweite (oft hn+1

2), typischerweise, bis

die gewünschte Toleranz erreicht istAlgorithmus:function [tsave,ysave] = _

adaptiveEinfach(f,y0,t0,T,h0,absTol)y = y0;t = t0;h = h0;% Paramatershmin = 1e-6;mu = 2;rho = 0.8;ysave = y0;tsave = t0;while t < T

k1 = f(t,y);k2 = f(t + h,y + h*k1);ytemp_Euler = y + h*k1;ytemp_Trapez = y + h/2*(k1+k2);if abs(ytemp_Trapez - ytemp_Euler) < absTol

y = ytemp_Trapez;t = t + h;% update h wie in Algorithmus 3.2 h = max(hmin, min(mu*h,rho*h*sqrt _ (absTol/norm(ytemp_Trapez - ytemp_Euler))));h = min(h, T - t);tsave = [tsave, t];ysave = [ysave, y];

elseh = h/2;if h < hmin

fprintf('h too small\n')return

endend

end

5.1 Schätzung des lokalen FehlersA posteriori Fehlerabschätzung, der im Schritt t n→t n+1 gemacht wurde:Berechne zuerst den Schritt mit der vorgegebenen Schrittweite. Danach

1. Halbierung der Schrittweite und Berechnung des Schritts nochmals mittels den Teilintervallen

2. Berechne den Schritt mit einem zweiten Verfahren höherer Ordnung

Verwende jeweils die Differenz der beiden Resultate zur Fehlerschätzung.

5.2 Toleranz des nächsten SchrittsZwei Strategien:

1. Versuche zu erreichen, dass der lokale Fehler proportional zu der Länge des Teilintervalls ist (d.h. für hk gross darf der Fehler auf

dem Intervall [ t k ,t k+1 ] größer sein als für hk klein)

z.B. ¿ lk+1=hk+1

T−t 0TOL

2. Versuche zu erreichen, dass der Fehler auf jedem Teilintervall gleich gross ist

Sonst auch eine Variante: Vorgeben von relativen und absoluten Toleranzen:

¿ lk+1=max ( reltol⋅‖yk‖,abstol )

5.3 Optimale SchrittweiteSeien yn+1

l ( yn ) und yn+ 1h ( yn ) die Werte zum Zeitpunkt t n+1 mit Hilfe

eines Verfahren niedriger Ordnung und eines höherer Ordnung.

Dann ist EST n+1=|yn+1h ( yn )− yn+1

l ( yn )| der geschätzte Fehler im

Schritt.Mache nun eine Probeschrittweite H und berechne basierend darauf die

Schrittweite hn+1¿

, die ab (t n , yn) eigentlich optimal wäre.Der Fehler der beiden Verfahren ist

e l= yn+1l ( yn )− y (t n; y n)≈c lH

p+1 und eh=…≈ chHp+2,

wenn angenommen wird, dass diese Ordnung p und p+1 haben (Konsistenzordnung ist eine Ordnung höher).Es kann gezeigt werden, dass somit

EST n+1∧≈|chHp+2−c lH

p+1|≈c lHp+1

¿⇔cl=ESTn+1

H p+1

und damit

c ⋅ (hn+1¿ )p+1≈¿

ln+1 ¿⇔hn+1¿ ≈ ( ¿ln+1

c )1

p+1 ≈c≔c l

H ⋅( ¿ ln+1

EST n+1)

1p+1 ¿

5.4 eingebettete RK-VerfahrenFalls zur Schätzung des Fehlers in einem Schritt zwei verschiedene RK-Verfahren mit unterschiedlicher Ordnung verwendet werden, ist man daran interessiert, dass so wenige Funktionsauswertungen wie möglich gemacht werden müssen. Eingebettete RK-Verfahren basieren darauf, dass die Stützstellen des Verfahrens mit tiefer Ordnung auch für jenes hoher Ordnung verwendet werden können.

Butcher –Tableau:

c Ac1

⋮cs

a11 … … a1 s

⋮ ⋱ ¿ ¿…¿…¿ ass¿

b lT b1

(l) ¿bs(l ) ¿

bhT b1

(h) ¿ bs(h )

ode45Basiert auf eingebetteten RK-Verfahren vierter und fünfter Ordnung.

5.5 Vergleich zu uniformer SchrittweiteI.d.R. genauer bei demselben numerischen Aufwand.Kann fehlschlagen, wenn grosse Anfangsfehler gemacht werden.



6 Numerische StabilitätFehler explodieren mit manchen Lösern: „Lösung ist nicht stabil“

6.1 Modellproblem und StabilitätsanalyseDAS Modellproblem für die Stabilitätsanalyse ist

y ( t )=λy ( t ) , y (t 0 )= y0 (¿ )

Die zugehörige analytische Lösung ist

y ( t )= y0 eλ ( t−t 0)

λ bestimmt das Verhalten von y ( t ) für t→∞.

Was passiert nun also für t→∞ ?

| y (t )∨¿|y0|⋅ exp( ℜ {λ }(t−t 0 ))={ ¿→0 ,∧¿ ℜ {λ }<0¿≡|y0|,∧¿ℜ {λ }=0¿→∞,∧¿ℜ {λ }>0

Wir wollen also ℜ {λ }<0.

Def. absolute Stabilität:Eine Einschrittmethode heisst absolut stabil für ein λh, wenn sie,

angewendet auf (*) für ℜ {λ }<0 beschränkte Lösungen erzeugt:

¿n≥0

¿ yn∨¿<∞¿

6.1.1 Stabilitätsfunktion / VerstärkungsfaktorFür die expl. Mittelpunktregel gilt gemäss Definition des Verfahrens:

{ k 1∧¿ f (t n , yn )

k2∧¿ f (t n+ h2 , yn+

h2 k1)

yn+1∧¿ yn+hk2

Einsetzen des Modellproblems, Faktoren zusammenfassen und nach g ( λh ) auflösen liefert:

yn+1∧¿ yn+hf ( yn+ h2f ( yn ))= yn+hf ( y n+

h2λ yn)

¿= yn+hλ( yn+h2λ yn)= yn+hλ yn+

h2 λ2

2yn= yn(1+hλ+ h2 λ2

2 )⏟g

⇔g (hλ )=1+hλ+ h2 λ2

2

⇒ Das Verfahren ist absolut stabil für |g (hλ )|<1

6.1.2 Stabilitätsgebiet / StabilitätsintervallStabilitätsgebiet:

SG≔ {z=hλ∈C :|g (z )|<1 }Stabilitätsintervall:

SI≔ {z=hλ∈R :|g ( z )|<1}Achtung! Es kann sein, dass g ( z ) nie z.B. negativ wird, dafür aber bei

mehreren Stellen ¿1.

6.1.2.1 explizite Verfahreng (hλ ) ist ein Polynom

Verfahren Stabilitätsfunktion Stabilitätsintervallexpl. Euler g ( z )=1+z (−2,0 )expl. Trapezregel

g ( z )=1+z+ z2

2

(−2,0 )

expl. Mittelpunktsregelg ( z )=1+z+ z2

2

(−2,0 )

klassisches RK4g ( z )=1+z+ z2

2+ z3

6+ z4

24

(−2.78,0 )

⇒ Handlungsrichtlinien für die Wahl von h:Das SG ist beschränkt, weshalb ein kleines h bei grossem λ verwendet werden muss.

6.1.2.2 Implizite Verfahren:g (hλ ) ist eine rationale Funktion

Verfahren Stabilitätsfunktion Stabilitätsgebiet

impl. Euler g ( z )= 11−z

impl. Trapezregel g ( z )=1+1

2z

1−12z

impl. Mittelpunktsregel g ( z )=1+1

2z

1−12z

Bei den in der Vorlesung vorgestellten Verfahren lässt sich feststellen, dass deren SG die ganze komplexe linke Halbebene abdeckt:

A-Stabilität: {z∈C : ℜ {z }<0 }⊂ SG



Für den impl. Euler gilt sogar noch g (−∞ )=0, was der Definition vonL-Stabilität entspricht.

⇒ Handlungsrichtlinien für die Wahl von h:h kann bei A-Stabilität beliebig gross gewählt werden (allerdings leidet die Genauigkeit).

6.1.2.3 Semi-Implizite VerfahrenDie Stabilitätsfunktionen der Semi-Impliziten Verfahren sind alle gleich wie jene der impliziten Verfahren weil das Newton-Verfahren in einem Schritt eine lineare Funktion (was das Modellproblem y ( t )=λy ( t ) ja auch ist) exakt berechnet.



6.2 BDF-Verfahrenlineare Mehrschrittmethode:für Berechnung von yn+ 1 werden auch yn , yn−1 ,… verwendet⇒ man muss keine Zwischenschritte auswerten

Nachteil:Schrittweitensteuerung ist sehr kompliziert

Vorteil:

A (α )-Stabilität: SG enthält einen Sektor {– α ≤π−arg (h )≤α }A (0 )-Stabilität: ganze negative reelle Achse ist im SG enthalten

6.3 Verwendung explizit / implizit / BDF

- implizite Verfahren haben bessere Stabilitätseigenschaften- implizite Verfahren sind teurer als explizite Verfahren- explizite Verfahren brauchen viel kleinere Schrittweiten als

implizite Verfahren für Stabilität - BDF: gut für reelle, negative Eigenwerte, da mindestens A (0 )-

stabil- ℜ { λ }>0 ⇒ explizite Verfahren verwenden

6.4 Systeme von ODE’sSei das folgende AWP gegeben:

y ( t )=A y ( t ) , y (0 )= y0 , A∈Cn×n

Vorausgesetzt, die Matrix A ist diagonalisierbar, lässt sich das ODE-System diagonalisieren. Es zerfällt in n entkoppelte ODE’s der Form

w i (t )=λ iwi (t )auf die jeweils einzeln die Stabilitätsanalyse gemacht werden kann.

6.5 Nichtlineare ProblemeUm eine Aussage zur (lokalen) Stabilität eines ESV machen zu können, muss das zugrunde liegende ODE-System linearisiert werden.Ein ESV mit dem SG ist dann numerisch stabil für ein gegebenes AWP, wenn die Schrittweiten hk so gewählt sind, dass für alle Eigenwerte λ i ( t ) der

Jacobi-Matrix DF (t , y (t )) (bzgl. des zweiten Arguments y (t )) mit

ℜ{λi (t ) }<0 gilt:

hk λ i (t k )∈SGDiese Hypothese ist mathematisch nicht rigoros, genügt aber, um im Allgemeinen gute Resultate zu erhalten.

Bsp:

Für das AWP

y (t )∧¿−200 t y2 (t )⏟f ( t , y (t ))

y (0 )∧¿1 gilt y ( t )= 1

1+100 t 2

⇒ f y (t , y (t ) )=−400 ty (t )=(¿ ) −400 t1+100t 2

(*): zum Glück kennt man hier die Lösungstrajektorie exakt.

Ein Plot von f y (t ) zeigt, dass das Minimum bei −20 liegt.

Für RK4 muss gelten hλ>−2.78⇔−20h←2.78

und damit h< 2.7820

≈0.139

6.6 Steife Problemey ( t )=A y ( t ) , y (0 )= y0 , A∈Cn× n

Betrachte die Eigenwerte von A:

κ=maxℜ {λi }<0

|ℜ {λ i }|min

ℜ {λi }<0|ℜ {λ i }|

≫1

Steif heisst für ein κ ab ca. 20.

Demnach:λ1=−2 , λ2=−100 ist steif

λ1=−200 , λ2=−100 ist nicht steif

Bei steifen Problemen sollten implizite Löser verwendet werden.



7 Appendix7.1 Jacobi-Matrixallgemein:

DF ( x )=(∇ f 1

T

∇ f 2T

⋮∇ f m

T )=(∂x1 f 1 ( x ) ∂ x2 f 1 ( x ) … ∂ xn f 1 ( x )∂x1 f 2 ( x ) ⋱ ¿ ¿

⋮ ¿⋱¿⋮ ¿∂x1 f m ( x )¿…¿…¿∂xn f 1 ( x )¿)2D:

DF ( x )=(∂ x1 f 1 ( x ) ∂ x2 f 1 ( x )∂x1 f 2 ( x ) ∂ x2 f 2 ( x ))

1D:

DF ( x )=∂x f (x )

7.2 Taylorallgemein:

j x0

n f (x )=∑k=0

n f ( k ) (x0 )k ! (x−x0 )k+O (|x−x0|

n+1 )Taylorentwicklungen:

y (t+h )= y ( t )+h⋅ y ' ( t )+ 12h2 ⋅ y ' ' ( t )+ 1

6h3⋅ y (3) (t )+O (h4 )

Berechnung bei Konsistenzordnung: y ( t ) und yk +1 sind Funktionen von h,

wenn bei y0 und k=0 begonnen wird. Entwickeln um h=0.Im Fall der Trapezregel:

y1 (h )= y0+h2¿

f (g (h ) ) ≈Taylor

∧f (g ( 0 ) )+h f ' (g (0 ) ) g' (0 )+ h2

2∂h [ f ' (g (0 ) )g ' (0 ) ]+O (h3 )

¿∧f ( y0 )+h f ' ( y 0 ) f ( y0 )+ h2

2f ' ' (g (0 ) ) f (g (0 ) ) f ( y0 )+O (h3 )

¿∧f ( y0 )+h f ' ( y0 ) f ( y0 )+ h2

2f ' ' ( y 0 ) f 2 ( y0 )+O (h3 )

mit

g (h )≔ y0+h f ( y0 ) , g (0 )=g0= y0 , g' (h )=g ' (0 )=f ( y0 )

7.3 Matlabfunction [retval1 retval2] = foo(argument)

for ii=0:99 %includes 0 and 99% Code here

endif expression

% Code hereendend

function expEulerConvf = @(t,y) y;t0 = 0;T = 2;y0 = 1;[~, exactsol] = ode45(f,[t0 T],y0);N = 2.^(0:10);err = [N zeros(1,length(N))];for ii = 1:length(N)

err(2,ii) = abs(expEuler(f,t0,T,y0,N(ii)) - ex-actsol(end));endloglog(err(1,:),’-*’,err(2,:),’-*’);legend(‚p = polyfit(log(err(1,:)),log(err(2,:)),1);fprintf('Ordnung der Quadratur ist %i\n', p(1));end

7.4 Finite DifferenzenTo do

7.5 WärmeleitungsgleichungTo do



7.6 Anmerkungen zur PrüfungAnforderungen an Zusammenfassung:5 Blätter (10 A4-Seiten) handgeschrieben. (Meine Abschrift dieses Dokuments war 7 Seiten lang)


Zusammenfassung Numerische Methoden http://blogs.ethz.ch/ricklis

1 Bemerkungen (V3.1)Diese Zusammenfassung ist aus meiner Mitschrift des PVKs Numerische Methoden von Ann-Katrin Tamm im Juni 2015 entstanden. Vielen herzlichen Dank für deine Arbeit, Ann-Katrin!

NotationFett gedruckte Symbole sind Vektoren, evtl. auch Matrizen.Gelbe Markierungen bezeichnen i.d.R. Abschnitte, welche eine Überarbeitung nötig haben, weitere Infos von den markierten Quellen nötig hätten oder unklar sind.

DisclaimerMeine Formelsammlungen entstehen und wachsen meist über eine längere Zeit. Es besteht immer ein gewisses Risiko, dass sich einige Fehler über zig Iterationen versteckt gehalten haben. Ich freue mich deshalb über jegliche (Fehler-) Verbesserungen, Anmerkungen, Lob, Dank oder auch Kritik.Meine Zusammenfassungen werden fortlaufend korrigiert und aktualisiert veröffentlicht.

Weiterverarbeitung:Weil ich es nicht ausstehen kann, dass ständig das Rad neu erfunden werden muss, habe ich das Originaldokument mit veröffentlicht mit der Einladung, sich hier für die eigene Formelsammlung zu bedienen. Ihr könnt diese Zusammenfassung also gerne weiterverarbeiten und / oder auch in überarbeiteter Form veröffentlichen.Haltet jedoch die Herkunft der kopierten/übernommenen Teile so gut wie möglich nachvollziehbar, falls ihr weiter veröffentlicht.Obiges gilt auch für alle anderen Formelsammlungen von mir, welche diesen Bemerkungstext (noch) nicht enthalten.

Quellenangaben:Aus Platz- und Zeitgründen (blame the 'Prüfungsstress') fehlen natürlich praktisch jegliche Quellenangaben (worüber ich auch schon ab und zu fluchen musste). Ich versuche jedoch in diesem Abschnitt die Arbeiten zu referenzieren, von denen ich wissentlich kopiert habe:

Wesentliche Bestandteile:Hauptbestandteil PVK von Ann-Katrin Tamm, Juni 2015Ergänzungen aus dem Skript von Dr. Sandra May

Revisionsverlauf:1.0 Sept 2015 erste Veröffentlichung Stefan Rickli1.1 28.09.2015 Notation vereinheitlicht Stefan Rickli

To Do: Finite Differenzen Wärmeleitungsgleichung

2 Ressourcen zu „Word und Formeleditor“Es gibt ein paar Ressourcen, welche mir sehr geholfen haben, den Formeleditor in Word zu meistern:

- Microsoft Word formula editor: https://support.office.com/en-us/article/Linear-format-equations-and-Math-AutoCor-rect-in-Word-2e00618d-b1fd-49d8-8cb4-8d17f25754f8?ui=en-US&rs=en-US&ad=US

- Unicode Nearly Plain-Text Encoding of Mathematics: http://www.unicode.org/notes/tn28/ o das ist der Standard, an den sich der Editor (fast vollständig) hält. Ist sehr gut beschrieben und

dokumentiert. Ausnahme sind Umrahmungen, welche nicht die ganze Funktionalität erhalten haben.- Die Zeichenübersicht des Editors selber: wenn man mit der Maus über ein Zeichen fährt, zeigt es einem den Tastatur-

Shortcut an, den man eingeben kann

Nice to know:- Alt + Shift + 0 erstellt eine neue Formel

(durch einen Bug in Office 2016 muss dieser Shortcut neu manuell definiert werden, Stand Mai 2016)- der Leerschlag ist euer Freund! Er veranlasst den Formeleditor, die Syntax bis zum aktuellen Punkt zu überprüfen und

das Zeug fixfertig bis zu dem Punkt, wo ihr seid, darzustellen (ausser es gibt noch offene Klammern).o verhält sich der Editor mal komisch, liegt es zu 95% daran, dass etwas in der Syntax nicht stimmt. Hier hilft

ab und an mal, sich die Formel im linearen Modus anzuschauen, wo alles bis auf Sonderzeichen wieder auseinander genommen wird. Das einzige, mit dem der Editor Mühe hat, sind grosse Eq-Arrays und mehrzeiliges Zeug.

- Wenn die Formel auf einer eigenen Zeile steht, veranlasst ein Leerschlag ausserhalb nach der Formel (also AUSSERHALB des Formeleditorfelds) den Editor, die Formel im Inline-Modus darzustellen (Formel braucht weniger Platz)

o siehe Tabellen in dieser ZF, dort habe ich das konsequent benutzt. Löscht mal das Leerzeichen nach einer Formel, das ein Integral enthält

- benutzt die Backslash (\) Befehle! Wenn man sich die beiden Dokumente oben ausdruckt und zur Referenz hält, geht es nicht lange, bis man alles mit der Tastatur machen kann und nie absetzen muss, um was mit der Maus zu machen

- \ensp und \emsp können helfen, um grössere, gewollte Abstände zu realisieren- \\eqarray ordnet mit jedem & einmal links und dann wieder rechts an- bastelt euch eure eigenen Shortcuts

o zum Beispiel \La für ⇐ \Ra für ⇒ \Lra für ⇔ oder \to für ein → oder eine leere 4x4 Matrix als \4x4 mit (■(&&&@&&&@&&&@&&&)) und einem Leerschlag

am Ende (damit der Ausdruck gleich aufgebaut wird) etc

o dazu einfach das entsprechende Zeichen in die Zwischenablage kopieren und im Formeleditor unter „Formeloptionen“ (in den Tools als kleiner Pfeil unten rechts zu finden) und „Math. Autokorrektur“ einfügen und den entsprechenden Backslashbefehl definieren

o manchmal ist es sinnvoll, noch eigene Funktionsnamen zu definieren, welche der Editor erkennen soll, wenn man sie häufig benutzt. Z.B. Imag()

ansonsten Leerschlag funktionsname\funcapply Leerschlag- die mathematische Autokorrektur ist manchmal auch ausserhalb des Formeleditors nützlich. Ich hab das in den

Optionen auch aktiviert- Wenn Word langsam wird wegen vielen anzuzeigenden Formeln, hilft

o 1. die Entwurfsansicht (anstatt Seitenlayout). Wenn man sich damit abfindet, dass dann ab und zu das Layout (noch) nicht dargestellt oder updated wird (nicht beirren lassen), kann man gut die kritischen Abschnitte bearbeiten. Achtung: Bilder werden NICHT angezeigt, sondern einfach mit einem Leerschlag repräsentiert!

o 2. reinzoomen, bis weniger Formeln sichtbar sind.

Zuletzt gespeichert: 03.01.2017, 19:42, Version 20 16 / 16

Persönliche Website – Services & Ressourcen | ETH …ricklis/blog/download... · Web...

Documents

Transcript of Persönliche Website – Services & Ressourcen | ETH …ricklis/blog/download... · Web...