Sekularna enacba, ki jo bomo resevali

f[x_] = 
      1 + 0.001/(1 - x) + 0.002/(2 - x) + 0.002/(3 - x) + 
        0.001/(4 - x)

Plot[f[x], {x, 0, 5}]

Pri zgornji sliki črta na sredini ni pri y=0, temveč pri y=1! Ničle so povsem blizu polov, zato se jih na zgornji sliki ne vidi.

Plot[f[x], {x, 0.9, 1.1}]

Točne reitve so naslednje. Očitno so ničle zelo blizu polom, zato Newtonova metoda ne bo delovala dobro. Morali bi imeti zelo dober začetni pribliek ampak tako dober začetni priblizek bi e skoraj pomenil, da poznamo tudi reitev. Mi pa bi radi imeli metodo, ki bo dobro delovala tudi, če vzamemo za začetni pribliek kar sredino med poloma (za tak začetni pribliek nam ni potrebno nič vedeti)

N[NSolve[f[x] == 0], 10]

Testi za Newtonovo metodo:

Newton[x_ ] := x - f[x]/f'[x];

Test[x_] := 
      Module[{y = {x}},  Do[AppendTo[y, N[Newton[Last[y]], 8]],  {5}];  
        y]

Poglejmo, kaj se dogaja z zaporedjem pri Newtonovi metodi, če napacno izberemo zaporedni pribliek. Funkcija Test izpie prvih pet priblikov:

Test[1.1]

Test[1.2]

Test[1.003]

Test[1.002]

Test[1.001]

Test[1.0001]

Test[2.01]

Test[2.005]

Test[2.004]

Test[2.002]

Vidimo, da mora biti začetni pribliek res zelo točen, sicer nam zaradi skoraj vodoravne tangente vse skupaj nekam pobegne.

Nov nastavek

Nov nastavek je, da delamo s psi1 in psi2, predpostavimo pa, da ičemo ničlo na intervalu med 2 in 3 (sicer sta psi1 in psi2 drugačna)

psi1[x_] = 0.001/(1 - x) + 0.002/(2 - x)
psi2[x_] = 0.002/(3 - x) + 0.001/(4 - x)

h1[x_] = c1/(2 - x) + cc1
h2[x_] = c2/(3 - x) + cc2

res1 = 
      NSolve[{h1[2.5] == 
            psi1[2.5], h1'[2.5] == psi1'[2.5]}, {c1, 
          cc1}]

res2 = 
      NSolve[{h2[2.5] == 
            psi2[2.5], h2'[2.5] == psi2'[2.5]}, {c2, 
          cc2}]

h[x_] = 
      First[First[1 + h1[x] + h2[x]  /.  res1  /.  res2]]

Plot[{f[x], h[x]}, {x, 2.0008, 2.02}]

N[NSolve[h[x] == 0], 10]

N[NSolve[f[x] == 0], 10]

Plot[f[x] - h[x], {x, 2, 3}]

Reitve sistema za b1, c1, b2 in c2 v tem primeru so

c1[x_] = psi1'[x] (2 - x)^2
cc1[x_] = psi1[x] - psi1'[x] (2 - x)
c2[x_] = psi2'[x] (3 - x)^2
\(cc2[x_] = 
      psi2[x] - \(\(\(psi2'\)[x]\) \((3 - x)\)\(
\)\(		\)\)\)

Za dani x izračunamo najprej c1,cc1,c2,cc2, nato pa izračunamo ničli funkcije h(x)=1+cc1+cc2+c1/(d_k-x)+c2/(d_(k+1)-x). Ičemo torej ničli kvadratnega polinoma, kar lahko vedno izračunamo direktno. Za začetni pribliek bomo vzeli 2.4.

Za kontrolo:

c1[2.4]
cc1[2.4]
c2[2.4]
cc2[2.4]

H[x_] := 
      1 + b1[x] + b2[x] + c1[x]/(2 - z) + c2[x]/(3 - z)

Racionalna funkcija, ki se ujema z vrednostjo f v tocki 2.4, prav tako pa se ujemata tudi odvoda, je :

h[z_] = H[2.4]

Poičemo ničle te racionalne funkcije. V resnici to naredimo tako, da jo spravimo na skupni imenovalec in poičemo ničle tevca, ki je kvadratni polinom.

N[Solve[h[z] == 0], 10]

Ena izmed dveh ničel lei znotraj pravega intervala, druga pa ne. Pravi naslednji pribliek je tista, ki lei na pravem intervalu, v naem primeru med 2 in 3. Postopek ponovimo tako, da vzamemo to za naslednjo točko.

h[z_]=H[2.002077425]

N[Solve[h[z] == 0], 10]

V dveh korakih smo iz 2.5 prili do skoraj povsem točnega rezultata.