Newtonin menetelmä ja iterointi

Kirjallisuutta

1. [KRE] Kreyszig: Advanced Engineering Mathematics, 8. painos Part D Complex Analysis, s. 651 - >
2. [BF] Burden-Faires
3. [RA] Kivelä: Reaalimuuttujan analyysi
4. [HubW] Hubbard-West
5. [Schei] Scheinerman: Invitation to dynamical systems
6. [HAM] Apiola: ... Maple ... (Newton s. 120, s. 65 - 67)

1. Yhtälöiden ratkaiseminen iteroimalla

1. Bisektio, perustuu Bolzanon lauseeseen [RA] Lause 3.4.2 s. 38.: Jos jatkuva funktio saa välin päätepisteissä erimerkkiset arvot, sillä on nollakohta välillä. Tässä on Matlab-funktio (keväältä 2000).
2. Newtonin menetelmä
3. Kiintopisteiteraatio
4. Interpolaatiomenetelmät
5. Erilaiset hybridimenetelmät

(1)                f(x) = 0

Bisektio

INPUT f,a,b,tol,maxiter
1. i:=1
2. WHILE i <= maxiter DO
     c:=(a+b)/2
     IF (f(c)=0 OR (b-a)/2 < tol) RETURN(c) FI
     i:=i+1
     IF f(a)*f(c) > 0, a:=c ELSE b:=c  # Tässä on asian ydin.
   ENDDO (OD)
   PRINT("Toleranssia ei saavuteta ",maxiter," "iteraatiolla")
END

          (b-a)/2n        

• Virhearvio ei riipu lainkaan funktiosta f, joten hitaasti ja varmasti
• Jos välillä on useita 0-kohtia, bisektioalgoritmi suppenee jotain kohti.

Yleinen iterointi

Alkupiste x0 .
          xn+1 = g(xn), n= 0,1,2,...

         g(x) = x .

Maple-välineitä iterointiin

> for k from 0 to 5 do x[k+1]:=g(x[k]) od:
> seq(x[k],k=0..5);
x[0], g(x[0]), g(g(x[0])), g(g(g(x[0]))), g(g(g(g(x[0])))),
    g(g(g(g(g(x[0])))))

seq((g@@k)(x[0]),k=0..5);
                      (2)           (3)           (4)           (5)
     x[0], g(x[0]), (g   )(x[0]), (g   )(x[0]), (g   )(x[0]), (g   )(x[0])

Kolme käytännön neuvoa (vanhan ja viisaan):

1. Ennenkuin kokeilet rakentamaasi iteraatiota, talleta työarkki.
2. Tee ensimmäinen kokeilu pienellä n:n arvolla
3. Muista käyttää evalf-komentoa. Iteraatioita tehdään liukuluvuilla, ei symboleilla (muidenkin on saatava käyttää laskentaresursseja).

Iteraation havainnollistus, aikasarja ja "cobweb"

Toinen hauska tapa on ns. "cobweb"-piirros, jossa piirretään samaan koordinaatistoon funktion kuvaaja sekä suora y=x. Näiden leikkauspisteitä haetaan (nehän ovat juuri kiintopisteitä).

Lähdetään jostain x₀:sta lasketaan y₀=f(x₀). Tämä on seuraava x₁, joten graafisesti se saadaan menemällä vaakajanaa y=y₀ suoralle y=x. Sitten vain jatketaan. Tässä Cobweb-opetusta (Univ. of Montana)

Esim. KRE s. 839 Exa 1

      x2 - 3x +1 =0

    g(x)=(1/3)(x2 +1)

  > g:=x->(1/3)*(x^2+1);
                                           2
                            g := x -> 1/3 x  + 1/3

  > X[0]:=1.0:n:=10:for i to n do X[i]:=g(X[i-1]) od:
  > seq(X[k],k=0..10);
1.0, .6666666666, .4814814814, .4106081389, .3895330145, .3839119898,
    .3824628053, .3820925991, .3819982514, .3819742213, .3819681019

> plot([seq([k,X[k]],k=0..10)]);  

Cobweb Elegantti Maple-toteutus saadaan määrittelemällä grafiikka-arvoinen funktio porras

 porras:=x->plot([[x,x],[x,g(x)],[g(x),g(x)]]):

 with(plots):
 g:=x->(1/3)*(x^2+1);
 X[0]:=1.0:n:=10:for i to n do X[i]:=g(X[i-1]) od:
 gjaxkuva:=plot([x,g(x)],x=0 .. 1,color=[blue,black]):
 alkup:=plot([[X[0],X[0]]],style=point,symbol=circle,symbolsize=20,color=blue):
 display([alkup,gjaxkuva,seq(porras(X[i]),i=0..n)]);

> X[10],g(X[10]);
                           .3819681019, .3819665436

Kiintopistelauseita

1. havainto

Jos g on jatkuva ja jos tiedämme, että kp-iteraatio x_n+1 = g(x_n) suppenee, sanokaamme c:tä kohti, niin c on g:n kiintopiste. Tämä seuraa heti jatkuvuudesta:

   xn+1=g(xn),

Tyypillinen klassinen "peruskurssiteknikka" rekursiivisesti määritellyn jonon raja-arvon laskemiseen on tällainen:

1. Saatetaan jono kiintopiste-iteraation muotoon, eli identifioidaan g,
2. Osoitetaan (induktiolla), että jono on kasvava ja ylhäältä rajoitettu (tai pienenevä ja alh. raj.) ==> on olemassa raja-arvo.
3. Ratkaistaan kp-yhtälö g(c)=c (jos osataan tai jos mahdollista). Ratkaisu c on kysytty raja-arvo.

Otamme samantien hieman tarkemmin kuin KRE-kirjassa, koska luonnollisista oletuksista saamme suorastaan ilmaiseksi myös kiintopisteen olemassaolo- ja yksikäsitteisyysjohtopäätöksen. (KRE-kirjassa oletetaan kiintopisteen olemassaolo.)

Tässä Latex-lähdekoodi. Taidanpa antaa Latex-kehyksen, jolla itsekukin voi ajaa nämä kauniisti. Toisaalta tätähän lukee kuin avointa kirjaa.

{\bf Kiintopistelause }

Olkoon $g\in C[a,b]$ ja $g(x)\in [a,b] \  \  \forall x\in [a,b]$, ts. $g$ kuvaa välin
$[a,b]$ sille itselleen. Tällöin $g$:llä on kiintopiste välillä $[a,b]$.
Jos lisäksi oletetaan, että $g$ on derivoituva ja
$$|g'(x)| \leq k \ \  \forall x\in [a,b],$$
missä $k<1$, niin kiintopiste on yksikäsitteinen ja iteraatiojono 
$x_n=g(x_{n-1}), n=0,1,\ldots$ suppenee mielivaltaisella alkupisteen $x_0$ 
valinnalla, $x_0\in [a,b]. $

Tod: 

(1) Olemassolo hoituu 

{\em Bolzanon lauseella}

Jos $g(a)=a$ tai $g(b)=b$, on kiintopiste 
löytynyt.
Oletetaan, että näin ei ole, eli $g(a)\ne a$ ja $g(b)\ne b$. Oletuksemme mukaan
täytyy silloin olla $g(a)>a$ ja $g(b)0$ ja  $g(b)<0$,
joten 
{\em Bolzano} $\implies \exists p\in [a,b]$ siten, että $h(p)=0$, eli 
$f(p)=p$. 

(2) Yksikäsitteisyys seuraa suoraan väliarvolauseesta:

Jos olisi kaksi eri kiintopistettä $p\ne q$ välillä $[a,b]$, niin
$$p-q=f(p)-f(q)=f'(\xi)(p-q),$$ joten saataisiin ristiriitainen epäyhtälö

$$|p-q| < k |p-q| < |p-q| .$$

(3) Suppeneminen seuraa epäyhtälöketjusta:

$|x_n -p| = |g(x_{n-1}) -g(p)| = |g'(\xi)||x_{n-1}-p| \leq k|x_{n-1}-p|
\leq \ldots \leq k^n |x_0-p|.$

Koska $0 < k < 1$, niin $\lim_{n\to\infty}x_n = p .$

Newton Maplella

Vrt. myös [HAM] ss. 120-121

> restart:
> N:=x->evalf(x-f(x)/D(f)(x));
> Nsymb:=x->x-f(x)/D(f)(x);
Käyttöesimerkki 
> f:=x->cos(x)-x;
> Nsymb(x); 
> seq((N@@n)(1.0),n=0..5);
> Matrix([seq([n,(N@@n)(1.0)],n=0..5)]); # Kätevä tapa taulukoida.








alaind:  
yläind: 
R^n   Rn