Musterlösung Aufgabe 19

Aufgabe 19: (Implementierungsaufgabe: das Newton-Maehly-Verfahren, 12 Punkte) Implementiere das Maehly-Verfahren zur automatischen Berechnung sämtlicher Nullstellen eines reellen Polynoms p(x) = a_{0}x^{n} + ⋯ + a_{n−1}x + a_{n} ohne komplexen Nullstellen. Eingabedaten: die Koeffizienten a_{0}, …, a_{n}, Ausgabedaten: alle Nullstellen x_{1} ≥ x_{2} ≥ … ≥ x_{n}. Intern sollen Funktionsauswertungen von p(x) und p′(x) über das Horner-Schema erfolgen.

a) Teste das Programm mit den Polynomen

\begin{matrix}p_{1}(x) = (x − 7)(x − 6) ⋯ (x − 1) = x^{7} − 28 x^{6} + 322 x^{5} − 1960 x^{4} + 6769 x^{3} \\ −13132 x^{2} + 13068 x − 5040 , \end{matrix}

p_{2}(x) = (x−3)(x−2)^{3}(x−1)(x+1)^{2} = x^{7}−8x^{6}+20x^{5}−6x^{4}−41x^{3}+38x^{2}+20x−24 .

b) Bestimme alle Eigenwerte der Matrix A = ( \begin{matrix}4 & 0 & 0 & 1 \\ 0 & 3 & 0 & 0 \\ 0 & 0 & 2 & 0 \\ 1 & 0 & 0 & α \end{matrix} ) mit α = M/7 777 777,

M= deine (bzw. eine eurer) Matrikelnummer(n). Das charakteristische Polynom von A ist

p_{3}(λ) = det(A − λ 1I ) = λ^{4} − (α + 9)λ^{3} + (9α + 25)λ^{2} − (26α + 19)λ + (24α − 6).

Beachte, daß A als symmetrische Matrix nur reelle Eigenwerte hat, sodaß p_{3} keine komplexen Nullstellen haben kann.

c) Betrachte das (eigentlich nicht zulässige) Polynom p_{4}(x) = p_{2}(x) + 0.001. Das Programm müßte alle 5 reellen Wurzeln finden. Versuche, die Lage der beiden nicht gefundenen Nullstellen in der komplexen Ebene mit Hilfe der gefundenen reellen Wurzeln zu ermitteln. (Hinweis: für p(x) = x^{n} + a_{1}x^{n−1} + ⋯ + a_{n} = ∏(x − z_{i}) folgt ∑ z_{i} = −a_{1} und ∏ z_{i} = (−1)^{n}a_{n}.)

Musterlösung: Lösung in MapleV.4:

horner:=proc(n,a,x)
local i,b,c:
b:=a[0]: c:=a[0]: for i from 1 to n-1 do b:=b*x+a[i]: c:=c*x+b: od:
[b*x+a[n], c];       # b*x+a[n] = polynom(x) , c = polynom'(x)
end:
maehly:=proc(n,a,anzahl_der_Nullstellen,genauigkeit)
local x_start,j,x,p,px,deltax,k,xx,
   Summe_der_restlichen_Wurzeln,
   Produkt_der_restlichen_Wurzeln:
Summe_der_restlichen_Wurzeln:=-a[1];
Produkt_der_restlichen_Wurzeln:=(-1)^n*a[n];
x_start:= evalf(max( abs(a[n]/a[0]), seq( 1+abs(a[i]/a[0]), i=1..n-1 ) )):
j:=1:
while j<=anzahl_der_Nullstellen do
   x:=x_start:
   while true do
     horner(n,a,x): p:=op(1,"); px:=op(2,"");
     deltax:=p/(px-p*sum(1/(x-xx[k]),k=1..j-1)):
     x:=x-deltax:
     if abs(deltax)<genauigkeit then break: fi:
   od:
   xx[j]:=x:
   Summe_der_restlichen_Wurzeln:=Summe_der_restlichen_Wurzeln-xx[j]:
   Produkt_der_restlichen_Wurzeln:=Produkt_der_restlichen_Wurzeln/xx[j]:
   j:=j+1;
od:
print(xx):
#  Summe  der restlichen Wurzeln muss Null sein
# Produkt der restlichen Wurzeln muss  1   sein
print(`  Summe der restlichen Wurzeln: `,  Summe_der_restlichen_Wurzeln):
print(`Produkt der restlichen Wurzeln: `,Produkt_der_restlichen_Wurzeln):
if anzahl_der_Nullstellen<n then
   print( `restliche Wurzeln: `,
   solve({a*b=Produkt_der_restlichen_Wurzeln,
          a+b=Produkt_der_restlichen_Wurzeln},{a,b})
        );
fi;
end:
Digits:=16:
# Aufgabenteil a), p1
n:=7 : a[0]:=   1:a[1]:= -  28: a[2]:=  322: a[3]:=-1960:
       a[4]:=6769:a[5]:=-13132: a[6]:=13068: a[7]:=-5040:

maehly(n,a,n,10^(-10)):
             table([ 1 = 6.999999999998736
                     2 = 5.999999999998538
                     3 = 4.999999999998433
                     4 = 4.000000000005333
                     5 = 2.999999999999919
                     6 = 2.000000000000080
                     7 = 1.000000000000000 ])
           Summe der restlichen Wurzeln:     -.1042*10^(-11)
           Produkt der restlichen Wurzeln: ,  .9999999999993915
# Aufgabenteil a), p2
n:=7 : a[0]:=  1:a[1]:= -8: a[2]:= 20: a[3]:= -6:
       a[4]:=-41:a[5]:= 38: a[6]:= 20: a[7]:=-24:
maehly(n,a,n,10^(-10)):
             table([ 1 =  2.999999999999969
                     2 =  1.999983815405862
                     3 =  1.999989595457301
                     4 =  2.000016725091570
                     5 =  1.000000000000000
                     6 = - .9999999952949274
                     7 = -1.000000004700197 ])
             Summe der restlichen Wurzeln: ,    .98640404224*10^(-5)
             Produkt der restlichen Wurzeln: , 1.000004932120923
# Aufgabenteil b), p3
alpha:=1.0*3777777/7777777:
n:=4:  a[0]:=1:            a[1]:= -9-  alpha:
       a[2]:= 25+ 9*alpha: a[3]:=-19-26*alpha:  a[4]:= -6+24*alpha:
maehly(n,a,n,10^(-10)):
                 table([ 1 = 4.264626419208630
                         2 = 2.999999999999994
                         3 = 1.999999999999997
                         4 =  .2210878150770824 ])
                 Summe der restlichen Wurzeln: ,    .56*10^(-14)
                 Produkt der restlichen Wurzeln: , 1.000000000000002
# Anmerkung: die exakten Wurzeln sind
exakt:= 3,2, 1/2*alpha+2 + 1/2*(alpha^2-8*alpha+20)^(1/2),
             1/2*alpha+2 - 1/2*(alpha^2-8*alpha+20)^(1/2);
                 exakt := 3, 2, 4.264626419208628, .221087815077082
# Aufgabenteil c), p4
n:=7 : a[0]:=  1:a[1]:= -8: a[2]:= 20: a[3]:= -6:

a[4]:=-41:a[5]:= 38: a[6]:= 20: a[7]:=-24+0.001:
maehly(n,a,n-2,10^(-10)):
table([ 1 =  2.999968746092964
2 =  2.047608993596442
3 = -1.002147619485920
4 = - .9978442785596099
5 =   .9998750390405419 ])
Summe der restlichen Wurzeln: ,   3.952539119315582
Produkt der restlichen Wurzeln: , 3.907411853551821
restliche Wurzeln: , 1.953705926775911 + .3007407608431249 I,
1.953705926775911 - .3007407608431249 I