H04 – Newton Verfahren, Lösen nichtlinearer Gleichungen

 

Schreiben Sie ein m-File („u41.m“), in dem Sie ein Newton-Verfahren zur iterativen Lösung einer nichtlinearen Gleichung f\left( x \right) = 0 programmieren. Die Funktion, die das Newton-Verfahren enthält, soll folgende Struktur haben:

function [x, it] = newton(x_0, tol, maxstep, f)
%
% x : approximierte Nullstelle
% it : Anzahl der Iterationen
% x_0 : Startwert der Iteration
% tol: Toleranz
% f : Funktion, deren Nullstelle bestimmt werden soll

Der Rückgabewert muss dann \left| {f\left( x \right)} \right| < tol erfüllen. Die als Parameter übergebene Funktion f wird in einem eigenen m-File realisiert. Sie bekommt die Auswertungsstelle(n) als Eingabe und gibt Funktionswert(e) und Ableitung(en) zurück. Der Aufruf erfolgt über newton(x_0, tol, maxstep, @f). Es hilft auch der feval-Befehl

Testen Sie das Programm für
a) {p_1}\left( x \right) = -{x^4}-2{x^3}+\frac{7}{2}{x^2}-\frac{{45}}{{16}},\quad {x_0} = 2
b) {p_2}\left( x \right) = {x^3}+4{x^2}+5x+2,\quad {x_0} = 0

Verwenden Sie zur Auswertung der Polynome das Horner-Schema. Beobachten Sie das Konvergenzverhalten, indem Sie sich den absoluten Fehler nach 1, 2, …, 10 Iterationen ausgeben lassen. Was fällt dabei auf? Wie können Sie die Konvergenz für {p_2} verbessern?

Lösung

Umstellen der Polynome:

{p_1}\left( x \right) = -{x^4}-2{x^3}+\frac{7}{2}{x^2}-\frac{{45}}{{16}} = -\frac{{45}}{{16}}+x\left( {\frac{9}{2}+x\left( {\frac{7}{2}-x\left( {2+x} \right)} \right)} \right)

{p_2}\left( x \right) = {x^3}+4{x^2}+5x+2 = 2+x\left( {5+x\left( {4+x} \right)} \right)

Ableitungen:

p_1^\prime \left( x \right) = -4{x^3}-6{x^2}+7x+\frac{9}{2} = \frac{9}{2}+x\left( {7-x\left( {6+4x} \right)} \right)

p_2^\prime \left( x \right) = 3{x^2}+8x+5 = 5+x\left( {8+3x} \right)

Matlab-Code

u41.m

format long;
[x1, it1] = newton(2, eps, 20, @p1)

newton.m

function [x, it] = newton(x_0, tol, maxstep, f)
%
% x: approx. Nullstelle
% it: Anzahl Iterationen
% x_0: Startwert der Iteration
% tol: Toleranz
% f: Funktion, deren Nullstelle bestimmt werden soll

x1 = x_0;
it = 0;
dev = zeros(10, 1);

gamma = 1;  % mit Korrektur für p2

while (it < maxstep) && (abs(f(x1)) > tol)
    it = it + 1;
    [y, dy] = f(x1);
    x = x1 - gamma * y / dy;
    if (it <= 10)
        dev(it) =  x - 1.5;
        % dev(it) = x + 1;
    end
    x1 = x;
end

dev

p1.m

function [y, dy] = p1(x)
y = -45/16 + x * (9/2 + x * (7/2 - x * (2 + x)));
dy = 9/2 + x * (7 - x * (6 + 4 * x));

p2.m

function [y, dy] = p2(x)
y = 2 + x * (5 + x * (4 + x));
dy = 5 + x * (8 + 3 * x);

Test für das erste Polynom

Verwendete Parameter: {x_0} = 2,\quad tol = eps,\quad maxstep = 20

Ergebnis: x = 1.5,\quad it = 6

Entwicklung des absoluten Fehlers:

\begin{array}{*{20}{c}}  i &\vline & {Abweichung} \\ \hline  1 &\vline & {{\text{0}}{\text{.185000000000000}}} \\  2 &\vline & {{\text{0}}{\text{.037988903850692}}} \\  3 &\vline & {{\text{0}}{\text{.002099903592718}}} \\  4 &\vline & {{\text{0}}{\text{.000006947947043}}} \\  5 &\vline & {{\text{0}}{\text{.000000000076432}}} \\  6 &\vline & {{\text{0}}{\text{.000000000000000}}} \\  7 &\vline & {{\text{0}}{\text{.000000000000000}}} \\  8 &\vline & {{\text{0}}{\text{.000000000000000}}} \\  9 &\vline & {{\text{0}}{\text{.000000000000000}}} \\{10} &\vline & {{\text{0}}{\text{.000000000000000}}} \\  \end{array}

Beobachtung: Die Konvergenz ist quadratisch.

Test 1 für das zweite Polynom

Verwendete Parameter: {x_0} = 0,\quad tol = eps,\quad maxstep = 20

Ergebnis: x = -0.999998252920598,\quad it = 20

Entwicklung des absoluten Fehlers:

\begin{array}{*{20}{c}}  i &\vline & {Abweichung} \\ \hline  1 &\vline & {{\text{0}}{\text{.600000000000000}}} \\  2 &\vline & {{\text{0}}{\text{.347368421052632}}} \\  3 &\vline & {{\text{0}}{\text{.193516663631397}}} \\  4 &\vline & {{\text{0}}{\text{.104014284855781}}} \\  5 &\vline & {{\text{0}}{\text{.054346842021410}}} \\  6 &\vline & {{\text{0}}{\text{.027856158851745}}} \\  7 &\vline & {{\text{0}}{\text{.014114290149138}}} \\  8 &\vline & {{\text{0}}{\text{.007105915824404}}} \\  9 &\vline & {{\text{0}}{\text{.003565448288787}}} \\{10} &\vline & {{\text{0}}{\text{.001785885343132}}} \\  \end{array}

Beobachtung: Die Konvergenz ist nur linear. Das liegt daran, dass die Steigung der Funktion in der Umgebung der Nullstelle sehr klein ist:

graph-polynom-newton

Zur Verbesserung der Konvergenz kann man das Newton-Verfahren abwandeln:

{x_{k+1}} = \tilde \Phi \left( {{x_k}} \right) = {x_k}-\gamma \frac{{f\left( {{x_k}} \right)}}{{{f^\prime }\left( {{x_k}} \right)}},\quad k \in \mathbb{N}

Dann gilt für die Konvergenz:

\tilde \Phi \left( x \right) = x-\gamma \frac{{\left( {x-\tilde x} \right)g\left( x \right)}}{{mg\left( x \right)+\left( {x-\tilde x} \right){g^\prime }\left( x \right)}}

{{\tilde \Phi }^\prime }\left( x \right) = 1+\frac{{\gamma g\left( x \right){g^{\prime \prime }}\left( x \right)\left( {x-\tilde x} \right)+m{g^\prime }\left( x \right)+{g^\prime }\left( x \right)\left( {x-\tilde x} \right)}}{{{{\left( {{g^\prime }\left( x \right)\left( {x-\tilde x} \right)+mg\left( x \right)} \right)}^2}}}

-\frac{{\gamma {g^\prime }\left( x \right)\left( {x-\tilde x} \right)}}{{{g^\prime }\left( x \right)\left( {x-\tilde x} \right)+mg\left( x \right)}}-\frac{{\gamma g\left( x \right)}}{{{g^\prime }\left( x \right)\left( {x-\tilde x} \right)+mg\left( x \right)}}

{{\tilde \Phi }^\prime }\left( {\tilde x} \right) = 1+\frac{{\gamma g\left( {\tilde x} \right)\left[ {{g^{\prime \prime }}\left( {\tilde x} \right)\left( {\tilde x-\tilde x} \right)+m{g^\prime }\left( {\tilde x} \right)+{g^\prime }\left( {\tilde x} \right)} \right]\left( {\tilde x-\tilde x} \right)}}{{{{\left( {{g^\prime }\left( {\tilde x} \right)\left( {\tilde x-\tilde x} \right)+mg\left( {\tilde x} \right)} \right)}^2}}}

\frac{{\gamma {g^\prime }\left( {\tilde x} \right)\left( {\tilde x-\tilde x} \right)}}{{{g^\prime }\left( {\tilde x} \right)\left( {\tilde x-\tilde x} \right)+mg\left( {\tilde x} \right)}}-\frac{{\gamma g\left( {\tilde x} \right)}}{{{g^\prime }\left( {\tilde x} \right)\left( {\tilde x-\tilde x} \right)+mg\left( {\tilde x} \right)}}

= 1-\frac{{\gamma g\left( {\tilde x} \right)}}{{mg\left( {\tilde x} \right)}} = 1-\frac{\gamma }{m}

Dabei ist m die Vielfachheit der Nullstelle.

Faktorisierung des gegebenen Polynoms:

{x^3}+4{x^2}+5x+2:\left( {x+1} \right) = {x^2}+3x+2

{x^2}+3x+2:\left( {x+1} \right) = x+2

\quad \Rightarrow \quad {x^3}+4{x^2}+5x+2 = \left( {x+1} \right)\left( {x+1} \right)\left( {x+2} \right)

Die Nullstelle x = -1 ist also doppelt. Die Konvergenz ist am schnellsten, wenn gilt:

{\tilde \Phi ^\prime }\left( {\tilde x} \right) = 0\quad \Rightarrow \quad 1-\frac{\gamma }{m} = 0\quad \Rightarrow \quad \gamma = m = 2

Test 2 für das zweite Polynom

Verwendete Parameter: {x_0} = 0,\quad tol = eps,\quad maxstep = 20,\quad \gamma = 2

Ergebnis: x = -{\text{0}}{\text{.999999993312917}},\quad it = 4

Entwicklung des absoluten Fehlers:

\begin{array}{*{20}{c}}  i &\vline & {Abweichung} \\ \hline  1 &\vline & {{\text{0}}{\text{.200000000000000}}} \\  2 &\vline & {{\text{0}}{\text{.015384615384614}}} \\  3 &\vline & {{\text{0}}{\text{.000115673799881}}} \\  4 &\vline & {{\text{0}}{\text{.000000006687083}}} \\  5 &\vline & {{\text{0}}{\text{.000000000000000}}} \\  6 &\vline & {{\text{0}}{\text{.000000000000000}}} \\  7 &\vline & {{\text{0}}{\text{.000000000000000}}} \\  8 &\vline & {{\text{0}}{\text{.000000000000000}}} \\  9 &\vline & {{\text{0}}{\text{.000000000000000}}} \\{10} &\vline & {{\text{0}}{\text{.000000000000000}}} \\  \end{array}

Beobachtung: Die Konvergenz des abgewandelten Verfahrens ist quadratisch.