Sovellus: Newtonin menetelmä
Lauri Ruotsalainen, 2010 (perustuu sovellukseen "Root Finding Using Bisection" by William Stein)
Seuraavassa ohjelmassa Newtonin algoritmin etenemistä kuvataan piirtämällä funktion tangenttisuoran suuntainen jana tarkastelupisteestä x-akselille. Janan ja vaaka-akselin leikkauspiste määrää seuraavan tarkastelupisteen paikan.
Menetelmän toimintaa voidaan tarkastella myös tulostamalla algoritmin suoritukseen liittyvien muuttujien arvot jokaisen iteraatiokierroksen jälkeen. Jälkimmäinen ohjelma tulostaa HTML-taulukon, jossa on listattuna iteraatiokierroksen järjestysnumero n, piste x_n, funktion arvo f(x_n) tässä pisteessä ja tulon f(c-h)f(c+h) arvo. Listauksesta nähdään, kuinka Newtonin menetelmän suoritus pysähtyy ensimmäisessä pisteessä, jossa f(c-h)f(c+h)<0.
Kuvia:
def newtonin_menetelma(f, c, maxn, h):
f(x) = f
valipisteet = [c]
kierros = 1
while True:
c = c - f(c)/derivative(f(x))(x=c)
valipisteet.append(c)
if f(c-h)*f(c+h) < 0 or kierros == maxn:
break
kierros += 1
return c, valipisteet
@interact
def _(f = x^2 - 2,
c = 6,
hh = slider(-16, -1, 1, -3, label="Tarkkuus 2h"),
maxn = slider(0, 15, 1, 10, label="Kierroksia max"),
tarkasteluvali = input_box(default = (0,6), label="Tarkasteluväli")):
f(x) = f
h = 10^hh
c, valipisteet = newtonin_menetelma(f, float(c), maxn, h/2.0)
html("$\\text{Tarkkuus = }%s$" %float(h))
html("$\\text{c = }%s$" %c)
html("$\\text{f(c) = }%s" %f(c))
html("$\\text{Iteratioita = }%s" %len(valipisteet))
# Piirretään algoritmin toimintaa kuvaava grafiikka.
P = plot(f, x, tarkasteluvali, color="blue")
L = sum(line([(c, 0), (c, f(c))], color="green") for c in valipisteet[:-1])
for i in range(len(valipisteet) - 1):
L += line([(valipisteet[i], f(valipisteet[i])), (valipisteet[i+1], 0)], color="red")
show(P + L, xmin=tarkasteluvali[0], xmax=tarkasteluvali[1], ymin=P.ymin(), ymax=P.ymax())
Click to the left again to hide and once more to show the dynamic interactive window |
def newtonin_menetelma(f, c, maxn, h):
f(x) = f
valipisteet = [c]
kierros = 1
while True:
c = c - f(c)/derivative(f(x))(x=c)
valipisteet.append(c)
if f(c-h)*f(c+h) < 0 or kierros == maxn:
break
kierros += 1
return c, valipisteet
@interact
def _(f = x^2 - 2, c = 6, hh = slider(-16, -1, 1, -3, label="Tarkkuus 2h"),
maxn = slider(0, 15, 1, 10, label="Kierroksia max")):
f(x) = f
h = 10^hh
c, valipisteet = newtonin_menetelma(f, float(c), maxn, h/2.0)
html("$\\text{Tarkkuus: } 2h = %s$"%float(h))
html("$\\text{Iteratioita = }%s<br>"%len(valipisteet))
# HTML-taulukon tulostaminen
s = "<br><br><table border=1 cellpadding=5>"
s += "<tr><td>$n$</td><td>$x_n$</td><td>$f(x_n)$</td><td>$f(x_n-h)f(x_n+h)$</td></tr>"
for i, c in enumerate(valipisteet):
s += "<tr><td>$%s$</td><td>$%s$</td><td>$%s$</td><td>$%s$</td></tr>"%(i, c, f(c), f(c-h)*f(c+h))
s += "</table>"
html(s)
Click to the left again to hide and once more to show the dynamic interactive window |
|
|