Laguerrova metoda

Izračun konstant s pomočjo programa Mathematica

S pomočjo programa Mathematica lahko elegantno izračunamo konstante pri razvoju ostanka Laguerrove metode po prejšnjih ostankih. Za začetek najprej na hitro ponovimo Laguerrovo metodo. Denimo, da iščemo korene funkcije f(z), z_{r} pa naj bo tekoči približek za enega od korenov. Naslednji približek potem izračunamo takole:

S_{1} = \frac{f′(z_{r})}{f(z_{r})},
S_{2} = \frac{f^{′2}(z_{r}) − f(z_{r})f′′(z_{r})}{f^{2}(z_{r})},
z_{r+1} = z_{r} − \frac{n}{S_{1} ± \sqrt{(n − 1)(nS_{2} − S_{1}^{2})}}.

V imenovalcu izberemo vedno tisti predznak, ki da imenovalcu večjo absolutno vrednost. Zgornjo formulo lahko spremenimo v obliko

z_{r+1} = z_{r} − \frac{nf(z_{r})}{f′(z_{r}) ± \sqrt{(n − 1)g(z_{r})}},  (1)

kjer je

g(z_{r}) = (n − 1)f^{′2}(z_{r}) − nf(z_{r})f′′(z_{r}).  (2)

Člene, ki nastopajo v zgornji formuli, razvijemo v Taylorjevo vrsto okrog korena a, najprej pa seveda še definiramo, da je f(a) = 0.

   In[1]:= f[a]=0
   In[2]:= F[z ]=Series[f[z],{z,a,3}]
   In[3]:= F'[z]=D[F[z],z]
   In[4]:= F''[z]=D[F'[z],z]

Tako smo v F[z], F'[z] in F''[z] shranili razvoje funkcije f in njenih odvodov okrog korena a. To sedaj vstavimo v definicijo (2) funkcije g.

   In[5]:= g[z ]=(n-1) F'[z]^2-n F[z] F''[z]

Nadaljujemo z razvojem \sqrt{(n − 1)g(z_{r})} v Taylorjevo vrsto.

   In[6]:= Series[Sqrt[(n-1) g[z]],{z,a,2}]

Ker Mathematica avtomatično ne poenostavi izraze oblike Sqrt[x^2] v x, jo moramo z ukazom PowerExpand spodbuditi k temu dejanju.

   In[7]:= PowerExpand[%]

Izraz se da še poenostaviti, kar vedno dosežemo z ukazom Simplify.

   In[8]:= Simplify[%]

Sedaj lahko izračunamo celotni razvoj imenovalca f′(z_{r}) + \sqrt{(n − 1)g(z_{r})}. Za predznak + se odločimo zaradi tega, ker bomo pri majhni absolutni vrednosti člena z − a tako dobili večjo absolutno vrednost v imenovalcu.

   In[9]:= F'[z]+%
   In[10]:= Simplify[%]

Ko imamo razvoj imenovalca v vrsto, lahko izračunamo celoten izraz (1).

   In[11]:= z-n F[z]/%
   In[12]:= Simplify[%]

Iz dobljenega rezultata lahko potem preberemo

ϵ_{r+1} = Cϵ_{r}^{3} + O(ϵ_{r}^{4}),

kjer je

C = \frac{(n − 2)f^{′′2}(α)}{8(n − 1)f^{′2}(α)} − \frac{f′′′(α)}{6f′(α)}.