Volker Traxler, TGM Wien XX
Newton-Raphson

Mathematische Inhalte:

Numerische Verfahren, Iteration. Anwendung: Newton-Verfahren. Kurzzusammenfassung: An Hand von Graphiken sollen Faustregeln für die Wahl eines Startpunktes bzw. Bedingungen für die Konvergenz des Newton-Verfahrens plausibel erläutert werden. Lehrplanbezug: 3. Jahrgang HTL; 7. Klasse AHS. Zeitaufwand: Je nach "Besprechungstiefe" 20 Minuten bis 3 Stunden. Mediales Umfeld: PCs; Overhead-LCD-Display; Mathematica Version 2.2 für Windows. Anmerkungen: Mit Hilfe der Software MathReader Version 2.2 (die an Schüler und Lehrer weitergegeben werden darf) können Sie den Inhalt meines Mathematica-Files, ohne im Besitz von Mathematica selbst zu sein, am PC betrachten und die Animationen vorführen !!!!. Allerdings können Sie keine Eingaben tätigen oder Berechnungen vornehmen !!!!

Bitte verwechseln Sie nicht MathReader (stark eingeschränkte Oberfläche, keine Eingabemöglichkeit etc.) mit dem CAS Mathematica !!!!

Der Ausdruck auf Papier gibt nicht die übliche Struktur eines Mathematica-Files in Form einer Inhaltsübersicht wieder und zeigt nicht die Möglichkeit des "Blätterns" in einem "elektronischen" Buch.

Da in Winword Graphik-Animationen nicht möglich sind, erscheint statt einer Graphik, die animiert werden kann, jeweils das letzte Einzelbild.

Die im Beispiel angeführten Formeln (I) bis (VII) wurden mit dem Formeleditor von Winword geschrieben und sind nicht Bestandteil des Mathematica-Files.

Newton-Raphson

1. Einführung

Für die folgenden Ausführungen wird vorausgesetzt, daß die Formel des Newton-Verfahrens bekannt bzw. abgeleitet und vielleicht schon an Hand einiger Beispiele (ohne auf das Problem der Wahl eines Startwertes einzugehen) der Algorithmus demonstriert wurde. Weiters kann angenommen werden, daß durch geschickte Auswahl von Beispielen bei der Selbsttätigkeit der Schüler Probleme durch ihre Auswahl der Startwerte auftraten.
Ziel der folgenden Ausführungen ist es, an Hand von Graphiken Faustregeln für die Wahl eines Startpunktes bzw. Bedingungen für die Konvergenz des Newton-Verfahrens (siehe Kapitel Zusammenfassung) plausibel zu machen.

2. Diskussion

n Wahl der Funktion f(x)

Für die folgenden Ausführungen bestimmen wir eine Funktion f(x), sowie die
dazugehörige Newtonsche Näherungsformel g(x).

In[13]:=

f[x_]:= x^3-4x+1; g[x_]:= x-f[x]/f'[x]

Stellen wir den Graphen der Funktion f(x) dar und speichern diese Graphik mittels der Variablen gr1.

In[14]:=

gr1 = Plot[f[x], {x,-3,3}];

Bemerkung: Der Vorteil eines "elektronischen" Notebooks besteht u.a. darin, daß der Schüler eine andere Funktion bzw. andere Startwerte wählen und die folgenden Befehle analog mit seinen Eingaben aufrufen kann. Er schafft sich damit seine eigene persönliche Mitschrift!

n Extremwerte der Funktion f(x)

Bei der Suche nach der Ursache von Problemen bei der Wahl eines geeigneten Startwertes wird meistens bald herausgefunden, daß f '(x) ungleich Null sein muß (Wiederholung: Definitionsmenge !!!). Die geometrische Interpretation der 1.Ableitung mußte öfters von den Schülern eingefordert werden.

Es folgt die Berechnung jener x-Werte, für die f ' (x) = 0 gilt:

In[15]:=

Solve[f'[x] == 0,x]

Out[15]:=

n Die Funktion PlotNewton

Um einen Eindruck weiterer Problemstellen zu bekommen, betrachten wir einige Animationen.
Bemerkung: mittels Eingabe von ? Befehlsname erhält man Hilfe, wie z. B.

In[16]:=

? PlotNewton

PlotNewton[f[x],{x,xmin,xmax},InitialPoint->x0,Iterations->m] plots f[x]
together with m Newton approximations starting at the initial point x0.
The sequence of approximations is also printed. The default value of
Iterations is 3. Graph->False suppresses the graph.
TangentStyle->style1 applies style1 to the tangent lines and points.
PointStyle->PointSize[n] sets the point size. Together->False produces
a series of plots suitable for animation. Other options are passed to
Plot.

In[17]:=

PlotNewton[f[x], {x,-3,3}, InitialPoint->-3,

Together->False];

{-3., -2.391304348, -2.154963544, -2.115945495}

In[18]:=

PlotNewton[f[x], {x,-3,3}, InitialPoint->-2/Sqrt[3],

Together->False];


Infinity::indet:
Indeterminate expression 1 + ComplexInfinity + ComplexInfinity
ncountered.
{-1.154700538, ComplexInfinity, Indeterminate, Indeterminate}

Mathematica weist uns mit einer Flut von Fehlermeldungen auf den Umstand hin, daß bei der obigen Rechnung der Nenner f '(x) Null und daher nicht definiert ist.

In[19]:=

PlotNewton[f[x], {x,-3,3}, InitialPoint->-.5,

Together->False];

{-0.5, 0.3846153846, 0.2492000512, 0.2540969477}

In[20]:=

PlotNewton[f[x], {x,-3,3}, InitialPoint->1,

Together->False];

{1., -1., 3., 2.304347826}

Aus obiger Animation entnimmt man: in der Nähe eines Extremwertes der Funktion f(x) kann es zu "Bocksprüngen" kommen.

In[21]:=

PlotNewton[f[x], {x,-3,3}, InitialPoint->3,

Together->False];

{3., 2.304347826, 1.967489477, 1.869470471}

Fazit: die Bereiche von -unendlich bis ca. -2. bzw. von ca. 1.8 bis + unendlich scheinen problemlos zu sein.

n Die Funktion xngraph

Um nun einen besseren Überblick zu bekommen, tragen wir in einem Diagramm eine Auswahl möglicher Startwerte x0 auf der x-Achse und als y-Wert den dazugehörigen z.B. 10-mal iterierten Wert ein, so erhalten wir folgendes Bild:

Diese Graphik kann problemlos mit folgendem Befehl xngraph (Name willkürlich, kann geändert werden) erzeugt werden.

Anforderung von Hilfe zu diesem selbstgestrickten Befehl erfolgt mit :

In[22]:=

? xngraph

Zeichnet eine Menge von Punkten mit den Koordinaten (Startwert/n-ter

iterierter Wert)

xngraph[g, xn, xmin, xmax, schrittweite, optionen]

g ..... eine Funktion (z.B g, nicht g[x] !!!)

xn .... Anzahl der Iterationen

xmin... linke Begrenzung der Graphik

xmax... rechte Begrenzung der Graphik

schrittweite

optionen ... Alle Optionen von ListPlot sind erlaubt!

Die so erzeugte Graphik wird unter dem Variablennamen gr2 abgelegt:

In[23]:=

gr2 = xngraph[g,10,-3,3,1/60];

Schon diese Graphik gibt Anlaß zu Diskussionen. Trotz 10-maliger Iterationen "erreichen" einige Startwerte nicht eine der drei reellen Wurzeln.

n Die Funktion xigraph

Tragen wir für verschiedene Startwerte wiederholt folgende Zahlenpaare

(1, x0), (2, g(x0) ), (3, g(g(x0)) ) in einem Koordinatensystem auf und legen die entstandene Graphik unter dem Variablennamen gr3 ab.

In[24]:=

? xigraph

Zeichnet eine Menge von Polygonzuegen
mit den Punkten (1,x0), (2,x1), (3,x2) .......(n+1,xn)
Eingabe analog xngraph!!!!

In[25]:=

gr3 = xigraph[g,10,-1.5,1.5,1/20];

Wie man sieht, spielt die Wahl eines geeigneten Startwertes eine entscheidende Rolle. Ein unglücklich gewählter Startwert führt entweder zu keinem Ergebnis oder zumindest zu einem gewaltigen Rechenaufwand.

n Fixpunkte von g(x)

Wir haben vorher bei der Betrachtung von g(x) bemerkt, daß f '(x) ungleich Null sein muß.
Erinnern wir uns an das Ziel unserer Betrachtungen: wir wollen jene x-Werte von f(x) bestimmen, für welche f(x) = 0 gilt.
Setzen wir dies in g(x) ein, so erhalten wir g(x) = x. Dies bedeutet, daß der x-Wert gleich dem y-Wert ist. (Wiederholung: 1. Mediane !!!!!!). Wir erhalten also einen Fixpunkt.

n Die Funktion graphik

Um diesen Sachverhalt zu demonstrieren, benutzen wir den Befehl graphik

In[26]:=

? graphik

Zeichnet den Graph der Funktion f(x), g(x) und der ersten Mediane.
f(x) erscheint schwarz, g(x) erscheint gruen, 1. Mediane blau
graphik[f,,xmin,xmax,{Pole von g(x)}, {Pole von f(x)}]
f ..... Funktion f (z.B f, nicht f[x] !!!)
xmin... linke Begrenzung, xmax... rechte Begrenzung
Pole von g(x) und f(x) in Form einer Liste, also in { } setzen !!!

In[27]:=

gr4 = graphik[f,-3,3, {-2/Sqrt[3],2/Sqrt[3]}];


Didaktischer Hinweis:Erzeugen Sie jeweils eine Folie mit dem Graph f(x), g(x) und der 1. Mediane und legen Sie diese nacheinander übereinander auf und verbinden Sie mittels Overheadstift die Punkte f(x) = 0 und die darüber liegenden Punkte g(x) = x. Die Verbindungslinien sind in Graphik gr4 nicht eingezeichnet.

n f(x) /f '(x)

Da wir ja i.a. keine exakten Lösungen der Gleichung f(x) = 0 angeben können, müssen wir erreichen, daß g(x) ungefähr gleich x wird, d.h. der Term |f(x) / f '(x)| soll möglichst klein sein.

Lassen wir uns f(x)/f '(x) graphisch darstellen.

In[28]:=

gr5 = PlotJump[f[x]/f'[x], {x,-3,3},

Jumps->{-2/Sqrt[3], 2/Sqrt[3]},

PlotStyle->RGBColor[.3,0,0]];


Überlagern wir die Graphiken gr1, gr2 und gr5:

In[29]:=

gr6 = Show[gr1,gr2,gr5, Frame->True];

Man sieht, daß in einer Umgebung der jeweiligen Nullstelle der Wert von |f(x) / f '(x)| tatsächlich klein ist , in der Nähe der Pole von f(x) / f ´(x) jedoch ein "Sprung" zu einer anderen Nullstelle erfolgt oder ein Erreichen der Nullstelle nach 10 Iterationsschritten noch immer nicht erfolgt ist.

Bemerkung: f(x) ungefähr gleich Null reicht i.a. nicht aus!!!!!!

n Die Funktion phasen

Wie nähern wir uns g(x) = x?
Wiederholen wir: wir nehmen einen (hoffentlich vernünftigen) Startwert x0, bilden g(x0), g(g(x0)), .... u.s.w.
Diese (hoffentlich konvergierende) Zahlenfolge kann mit dem Mathematica-Befehl NestList erzeugt werden.

In[30]:=

? NestList

NestList[f, expr, n] gives a list of the results of applying f to expr 0
through n times.

In[31]:=

liste1 = NestList[g,3,3] // N

Out[31]:=

{3., 2.30435, 1.96749, 1.86947}

Vergleiche mit der Ausgabe von PlotNewton mit Startwert 3 weiter oben!
Wie weiter oben bereits angerissen, wird das Finden einer Nullstelle einer Funktion f(x) zurückgeführt auf das Finden eines Fixpunktes einer Funktion g(x). Geometrisch interpretiert ist dies der Schnittpunkt der ersten Mediane mit g(x) (siehe Graphik gr4).
Mit dieser, sogenannten Fixpunktiteration, bilden wir eine Folge von Punkten,

(x0, g(x0)), (g(x0), g(g(x0))) u.s.w.

Mit Startwert 3 erhalten wir für die oben definierte Funktion g(x) die Punktfolge

(3, 2.30435), (2.30435, 1.96749), (1.96749, 1.86947)

Versuchen wir diesen Sachverhalt für eine Reihe von verschiedenen Startwerten graphisch darzustellen:

In[32]:=

? phasen

Zeichnet das Phasendiagramm der Iterationsfunktion g(x)
phasen[g, xn, x0min, x0max, schrittweite, optionen]
g ..... Funktion g (z.B g, nicht g[x] eingeben !!!)
xn ..... Anzahl der Iterationen
x0min .. kleinster Startwert
x0max .. größter Startwert
schrittweite
optionen siehe Optionen von ListPlot:

In[33]:=

gr7 = phasen[g, 5, -3,3,1/60];

Diese Graphik kommt uns bekannt vor. Sie entspricht der Graphik der Funktion g(x) (siehe gr4). Dabei beobachtet man die Anhäufung der einzelnen Iterationspunkte bei den Fixpunkten der Funktion g(x). (Stichwort: Grenzwert !!!)

Stellen wir die Folge der Zahlenwerte (3, 2.30435), (2.30435, 1.96749), (1.96749, 1.86947) dar:

In[34]:=

gr8 = phasen[g,10,3,3.5,3, PlotStyle-> {RGBColor[0,.5,.5],

PointSize[.03]}];

 

n Die Funktion iteration

Gibt es eine graphische Möglichkeit, rasch diese Punktfolge zu zeichnen ?
Man kann die einzelnen Punkte mittels eines "Zick-Zack Bandes" zwischen dem Graphen von g(x) und der ersten Mediane verbinden. Dazu kann man den Befehl iteration benutzen:

In[35]:=

? iteration

Geometrische Darstellung der Iteration einer Funktion g als Zick-Zack-Band zwischen dem Graphen der Funktion g und der 1. Mediane!
iteration[g, x0, xi, xn, xmin, xmax, pole]
g .... zu iterierende Funktion
x0 .... Startwert
xi .... Zwischenwert, ab dem die Iterationen dargestellt werden sollen
xn .... Endwert
xmin .. untere Grenze des Zeichenbereichs
xmax .. obere Grenze des Zeichenbereichs
pole .. Liste der Pole der Funktion
Achtung: Pole in Form einer Liste eingeben, z.B.: {-7, 0.1, 3}

In[36]:=

gr9 = iteration[g,3,0,10,1.5,3];

Überlagerung der Graphiken gr8 und gr9 ergibt:

In[37]:=

gr10 = Show[gr8, gr9, PlotRange->All];

Daß die Iterationspunkte auf g(x) liegen, zeigt die Graphik gr11:

In[38]:=

gr11 = Show[gr4,gr7, PlotRange->{-7.5, 7.5}];

Wie sieht es nun mit der Konvergenz der Iteration für einzelne Startpunkte aus? Weiter oben (siehe PlotNewton) wurde bereits vermutet, daß die Bereiche von -unendlich bis ca. -2 bzw. von ca. 1,8 bis + unendlich problemlos zu sein scheinen.
Leicht können andere Startwerte ausprobiert werden.

n |g ' (x)| < 1

In[39]:=

gr12 = iteration[g,1.2,0,10,-1.5,8, {-2/Sqrt[3],2/Sqrt[3]}];


Wie man sieht, ist der rechte obere Zweig von g(x) eine konvexe Funktion (Stichwort: Kurvendiskussion). Offensichtlich gilt in einer Umgebung des Fixpunktes (1.86, 1.86), daß g'(x) < 1 ist. Genauere Untersuchungen zeigen, daß -salopp ausgedrückt- in einem Intervall um den Fixpunkt |g'(x)| < 1 gelten muß.

Definieren wir in Analogie zu f(x) und g(x) zwei Funktionen u(x) und v(x), um die Definitionen für f(x) und g(x) nicht zu überschreiben:

In[40]:=

v[x_] := x-u[x]/u'[x]

Die erste Ableitung ergibt:

In[41]:=

v'[x]

Out[41]:=
 

u[x] u''[x]
-----------
     2
u'[x]
Offensichtlich muß f(x) mindestens zweimal differenzierbar sein.

n f ''(x)

Untersuchen wir die Rolle der zweiten Ableitung von f(x):

In[42]:=

gr13 = Plot[Evaluate[{f[x], f''[x]}, {x,-3,3},

PlotStyle->{{},{RGBColor[0,0.5,.5]}}]];

In[43]:=

gr14 = Show[gr2,gr13];

Man erkennt: wenn sgn f(x) = sgn f '' (x) gilt, konvergiert die Iterationsfolge mit einem Startwert aus diesem Bereich sicher zu der diesem Bereich nächstgelegenen Nullstelle.

3. Zusammenfassung

Es ist also unbedingt notwendig, den ungefähren Bereich der Nullstelle graphisch zu bestimmen.
Für die Randpunkte a und b des Bereichs sollte gelten: f(a)* f(b) < 0
Die Funktion sollte in einer Umgebung der Nullstelle zweimal stetig differenzierbar sein
Die erste Ableitung darf für alle Werte aus dieser Umgebung nicht Null ergeben
Der Startwert muß hinreichend "nahe" bei der Nullstelle liegen
Praxis: sgn f(x) = sgn f''(x)

4. Beispiel

Für ein Bohrgestänge benötigt man ein Rohr mit einem Innendurchmesser von . Um der Belastung standzuhalten, soll das Rohr ein Drehmoment aufnehmen können. Die zulässige Torsionsspannung beträgt.

Welcher Außendurchmesser  ist erforderlich?

Es gilt: (siehe Buchliste: [1] Seite 233)

(I)  und (II)  und (III) ,

wo  das polare Flächenmoment und  das polare Widerstandsmoment darstellt. 

und  sind der jeweilige Innen- bzw. Außendurchmesser.

Aus (I) folgt

.

Die technische Anforderung an das polare Widerstandsmoment ergibt sich daher zu

(IV)     .

Aus (III) folgt

(V)     

(V) eingesetzt in (II) ergibt:

(VI)     

Mit den Bedingungen (IV) und (VI) ergibt sich die folgende Ungleichung

(VII)     

Umformung und Einsetzen des Wertes für den Innendurchmesser liefert die Ungleichung für den erforderlichen Außendurchmesser

(VIII)     

Anmerkung: Ein Teil der Schüler erkennt hier nicht die Notwendigkeit eines numerischen Verfahrens ("numerisch? Geht ja nicht, ist ja eine Ungleichung") und ein Teil hat bei der "händischen" Suche nach einen vernünftigen Startwert interessanterweise Probleme ("wie soll ich bei so großen Zahlen einen vernünftigen Startwert finden?").

 

In[44]:=

Plot[Pi da^4-Pi 50^4-16 20000 da, {da,40,100}];

Mit dem Innendurchmesser als Startwert sind alle in der Zusammenfassung angeführten Faustregeln erfüllt.

Wir erhalten:

In[45]:=

PlotNewton[ Pi da^4-Pi 50^4-16 20000 da, {da,50,61},

InitialPoint->50];

{50., 62.79185081, 59.52786986, 59.19978184}

Zur Kontrolle die Lösung mit einem in Mathematica eingebauten Befehl:

In[46]:=

NSolve[Pi da^4-Pi 50^4-16 20000 da == 0,da]

Out[47]:=

{{da -> -38.8937}, {da -> -10.1515 - 51.1032 I},

{da -> -10.1515 + 51.1032 I}, {da -> 59.1966}}

Das Rohr muß einen Außendurchmesser von 60 mm aufweisen.

Hinweis:
Folgende Befehle sind vom Autor selbst entwickelt worden und sind nicht im Standardpaket Mathematica enthalten: xngraph, xigraph, graphik, phasen, iteration.

Die Befehle PlotNewton und PlotJump wurden aus einem ebenfalls nicht im Standardpaket enthaltenen Package ("KnoxPackage", siehe Anweisungen auf Diskette) eingelesen.

Bei Aufruf des Programms durch Mathematica werden durch Anklicken der Ja-Button all diese Befehle (sofern KnoxPackage ordnungsgemäß installiert) automatisch eingelesen.

Buchliste:

[1] Schärf, Julius : "Mathematik 3 für HTL", 6. Auflage, R. Oldenburg Verlag Wien 1994,
ISBN 3-7029-0389-5