Newton-Raphson |
Mathematische Inhalte:
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.
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]Offensichtlich muß f(x) mindestens zweimal differenzierbar sein.
-----------
2u'[x]
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: