Grundlagen der funktionsorientierten Programmierung mit SCHEME/Das Newton-Verfahren
Aus ProgrammingWiki
Wiki-Quelltext dieser Seite.
<language>Scheme</language>
'''Autoren: Schneider, Tim; Berger, Veit'''
== Grundlagen ==
Zur näherungsweisen Bestimmung der Nullstelle einer stetigen Funktion <math>f</math> kann an einer beliebigen Stelle <math>x</math> die Tangente angelegt werden. Für den Anstieg <math>m</math> der Tangente gilt:
<blockquote><math>m = f'(x) = \frac{f(x)}{x-x_0}</math></blockquote>
wobei <math>x_0</math> die Nullstelle der Tangente ist. Diese Nullstelle ist ein erster Näherungswert für die Nullstelle der Funktion <math>f</math> und wird als Startwert für die erneute Anwendung des Verfahrens genutzt.
Die Nullstellensuche soll abgebrochen werden, wenn <math>|f(x)|</math> in einer hinreichend kleinen <math>\varepsilon</math>-Umgebung liegt.<br />
Damit kann die Nullstelle der Funktion <math>f</math> mit beliebiger Genauigkeit angenähert werden.
== Algorithmus ==
<ul>
<li>
Lege einen geeigneten Startwert <math>x</math> fest.
</li>
<li>
Falls <math>|f(x)|<\varepsilon</math> ist, dann gib <math>x</math> als Nullstelle zurück.
</li>
<li>
Anderenfalls berechne
<blockquote><math>x_0 = x - \frac{f(x)}{f'(x)}</math></blockquote>
und suche die Nullstelle mit dem neuen Startwert <math>x_0</math>.
</li>
</ul>
== Implementation ==
Zur Bestimmung des Anstiegs <math>m=f'(x)</math> nutzen wir die Prozedur höherer Ordnung <tt>ableitung</tt> (s. [[Grundlagen der funktionsorientierten Programmierung mit SCHEME/Prozeduren höherer Ordnung (II)|Prozeduren höherer Ordnung (II)]]).<br />
Natürlich müssen wir auch daran denken, dass die stetige Funktion <math>f</math> u. U. keine reelle Nullstelle besitzt.
<code id="53ea65952db7c" hidden>
(define sqr
(lambda (x)
(* x x)))
</code>
<code id="53ea65952dcde">
(define h 1e-10)
(define eps 1e-08)
(define maxanz 100)
(define ableitung
(lambda (fkt)
(lambda (x)
(/ (- (fkt (+ x h)) (fkt x)) h))))
(define nst
(lambda (fkt x)
(letrec ([suchen
(lambda (x anz)
(cond
[(< anz 1) #f]
[(< (abs (fkt x)) eps) x]
[else
(suchen (- x (/ (fkt x) ((ableitung fkt) x))) (- anz 1))]))])
(suchen x maxanz))))
</code>
<run id="53ea65952e06e">
(nst (lambda (x) (+ (expt x 3) (* -2 (expt x 2)) (* -11 x) 12)) 3)
</run>
== Nullstellensuche mit Schaubild ==
Wir wollen die Nullstellensuche am Schaubild der Funktion demonstrieren und damit das Ergebnis überprüfen. Dazu sammeln wir mit der Prozedur <tt>fkt-list</tt> die Lambda-Terme der Funktion und die aller Tangenten, die zur Nullstellensuche benötigt werden.
Abschließend wollen wir das numerische Ergebnis sinnvoll runden.
<code id="53ea65952de04" hidden>
(define nr 1)
(define genau 6)
(define breite 610)
(define hoehe 460)
(define rand 50)
(define canvas (eval (string->symbol (string-append "canvas" (number->string nr)))))
(define raster #f)
(define bez "x")
(define black '(0 0 0))
(define red '(255 0 0))
(define brown '(100 50 0))
(define orange '(255 100 0))
(define yellow '(255 200 0))
(define green '(0 255 0))
(define blue '(0 0 255))
(define pink '(200 50 200))
(define rgb-list
(lambda (anz)
(letrec
([rgb
(lambda (n ls)
(if (< n 1)
(cons '(255 0 0) ls)
(rgb (- n 1) (cons `(0 ,(round (- 255 (* n (/ 255 anz)))) ,(round (- 255 (* n (/ 255 anz))))) ls))))])
(rgb (- anz 1) ()))))
(define n-tes-mod
(lambda (ls anz)
(list-ref ls (modulo anz (length ls)))))
(define genaurunden
(lambda (x anz)
(let ([n 0] [erg 0])
(if (not (= x 0))
(set! n (inexact->exact (round (- anz (/ (log (abs x)) (log 10)) -0.5)))))
(set! erg (/ (round (* x (expt 10 n))) (expt 10 n)))
(if (integer? erg)
(inexact->exact erg)
(exact->inexact erg)))))
(define groessenordnung
(lambda (x)
(if (= x 0)
1
(let ([pot (/ (log (abs x)) (log 10))])
(genaurunden (expt 10 (floor pot)) genau)))))
(define minimum
(lambda (x y)
(if (< x y) x y)))
(define maximum
(lambda (x y)
(if (> x y) x y)))
(define ungleich-null
(lambda (x y)
(if (= x y 0) (+ y 1) y)))
(define schaubild
(lambda (fkt-list x1 x2 y1 y2 von bis koordcolor color-list)
(clear canvas)
(setwidth canvas 1)
(let* ([textgroesse 9]
[rand 50]
[xu (exact->inexact (minimum x1 x2))] [xo (exact->inexact (maximum x1 x2))] [xo (ungleich-null xu xo)]
[yu (exact->inexact (minimum y1 y2))] [yo (exact->inexact (maximum y1 y2))] [yo (ungleich-null yu yo)]
[spmin rand] [spmax (- breite rand)] [spend (+ spmax 15)]
[zemin (- hoehe rand)] [zemax rand] [zeend (- zemax 15)]
[s00 (+ (* (/ (- xu) (- xo xu)) (- spmax spmin)) spmin)] [s0 s00]
[z00 (+ (* (/ (- yu) (- yo yu)) (- zemax zemin)) zemin)] [z0 z00]
[spvon (+ (* (/ (- von xu) (- xo xu)) (- spmax spmin)) spmin)]
[spbis (+ (* (/ (- bis xu) (- xo xu)) (- spmax spmin)) spmin)]
[mx (- xo xu)] [teilerx 1] [stepx (/ (- spmax spmin) mx teilerx)] [ex mx] [faktorx 1]
[my (- yo yu)] [teilery 1] [stepy (/ (- zemin zemax) my teilery)] [ey my] [faktory 1])
(letrec
([anpassen
(lambda (m)
(if (> m 2)
(if (<= m 20)
m
(anpassen (/ m 10)))
(if (> m 2)
m
(anpassen (* m 10)))))]
[teilen
(lambda (m)
(cond
[(< m 3) 5]
[(< m 10) 2]
[else 1]))]
[x-raster
(lambda (sx zy)
(if (> zy (- zemax 1))
(begin
(if (and (< zy (+ zemin 1)) (> (abs (- zy z0)) 1)) (line canvas sx zy sx zy))
(x-raster sx (- zy 3)))))]
[x-teilen-positiv
(lambda (sx)
(if (< sx (+ spmax 1))
(begin
(if (and (> sx (- spmin 1)) (> (abs (- sx s0)) 1))
(begin
(line canvas sx z0 sx (- z0 3))
(if raster (x-raster sx zemin))))
(x-teilen-positiv (+ sx stepx)))))]
[x-teilen-negativ
(lambda (sx)
(if (> sx (- spmin 1))
(begin
(if (and (< sx (+ spmax 1)) (> (abs (- sx s0)) 1))
(begin
(line canvas sx z0 sx (- z0 3))
(if raster (x-raster sx zemin))))
(x-teilen-negativ (- sx stepx)))))]
[y-raster
(lambda (sx zy)
(if (< sx (+ spmax 1))
(begin
(if (and (> sx (- spmin 1)) (> (abs (- sx s0)) 1)) (line canvas sx zy sx zy))
(y-raster (+ sx 3) zy))))]
[y-teilen-positiv
(lambda (zy)
(if (> zy (- zemax 1))
(begin
(if (and (< zy (+ zemin 1)) (> (abs (- zy z0)) 1))
(begin
(line canvas s0 zy (+ s0 3) zy)
(if raster (y-raster spmin zy))))
(y-teilen-positiv (- zy stepy)))))]
[y-teilen-negativ
(lambda (zy)
(if (< zy (+ zemin 1))
(begin
(if (and (> zy (- zemax 1)) (> (abs (- zy z0)) 1))
(begin
(line canvas s0 zy (+ s0 3) zy)
(if raster (y-raster spmin zy))))
(y-teilen-negativ (+ zy stepy)))))]
[faktor
(lambda (len step)
(if (<= len step)
1
(+ 1 (faktor (- len step) step))))]
[textwidth
(lambda (str)
(* textgroesse (string-length str)))]
[textout-widthcenter
(lambda (sx zy str unten)
(if unten
(if (integer? (string->number str))
(text canvas (- sx (/ (textwidth str) 2))
(+ zy (/ textgroesse 2/3)) str textgroesse)
(text canvas (- sx (/ (textwidth str) 2) (- (/ textgroesse 2)))
(+ zy (/ textgroesse 2/3)) str textgroesse))
(if (integer? (string->number str))
(text canvas (- sx (/ (textwidth str) 2))
(- zy textgroesse) str textgroesse)
(text canvas (- sx (/ (textwidth str) 2) (- (/ textgroesse 2)))
(- zy textgroesse) str textgroesse))))]
[textout-heightcenter
(lambda (sx zy str links)
(if links
(if (integer? (string->number str))
(text canvas (- sx (textwidth str) 2)
(+ zy (/ textgroesse 2)) str textgroesse)
(text canvas (- sx (textwidth str) (- (/ textgroesse 2)))
(+ zy (/ textgroesse 2)) str textgroesse))
(if (integer? (string->number str))
(text canvas (+ sx textgroesse)
(+ zy (/ textgroesse 2)) str textgroesse)
(text canvas (+ sx textgroesse)
(+ zy (/ textgroesse 2)) str textgroesse))))]
[textout-bez-widthcenter
(lambda (sx zy str)
(text canvas (- sx (/ (textwidth str) 4))
(- zy textgroesse) str textgroesse))]
[textout-bez-heightcenter
(lambda (sx zy str)
(text canvas (+ sx (/ textgroesse 2))
(+ zy (/ textgroesse 2)) str textgroesse))]
[x-bezeichnen-positiv
(lambda (sx e unten)
(if (< sx (+ spmax 1))
(begin
(if (and (> sx (- spmin 1)) (> (abs (- sx s0)) 1))
(textout-widthcenter sx z0 (number->string (genaurunden e genau)) unten))
(x-bezeichnen-positiv (+ sx (* stepx faktorx)) (+ e (* ex faktorx)) unten))))]
[x-bezeichnen-negativ
(lambda (sx e unten)
(if (> sx (- spmin 1))
(begin
(if (and (< sx (+ spmax 1)) (> (abs (- sx s0)) 1))
(textout-widthcenter sx z0 (number->string (genaurunden e genau)) unten))
(x-bezeichnen-negativ (- sx (* stepx faktorx)) (- e (* ex faktorx)) unten))))]
[y-bezeichnen-positiv
(lambda (zy e links)
(if (> zy (- zemax 1))
(begin
(if (and (< zy (+ zemin 1)) (> (abs (- zy z0)) 1))
(textout-heightcenter s0 zy (number->string (genaurunden e genau)) links))
(y-bezeichnen-positiv (- zy (* stepy faktory)) (+ e (* ey faktory)) links))))]
[y-bezeichnen-negativ
(lambda (zy e links)
(if (< zy (+ zemin 1))
(begin
(if (and (> zy (- zemax 1)) (> (abs (- zy z0)) 1))
(textout-heightcenter s0 zy (number->string (genaurunden e genau)) links))
(y-bezeichnen-negativ (+ zy (* stepy faktory)) (- e (* ey faktory)) links))))]
[x-achse-bezeichnen
(lambda (str)
(textout-bez-heightcenter spend z0 str))]
[y-achse-bezeichnen
(lambda (str)
(let ([y-str (string-append "f(" str ")")])
(textout-bez-widthcenter s0 zeend y-str)))]
[graph-startpkt
(lambda (fkt pkt)
(let ([sx (car pkt)] [zy (cdr pkt)])
(if (< sx (+ spbis 1))
(with-failure-continuation
(lambda (error-record error-k) (graph-startpkt fkt (cons (+ sx 1) zy)))
(lambda ()
(let* ([x (+ (* (/ (- sx spmin) (- spmax spmin)) (- xo xu)) xu)] [y (fkt x)])
(if (real? y)
(begin
(set! zy (+ (* (/ (- y yu) (- yo yu)) (- zemax zemin)) zemin))
(if (> zy zemin) (set! zy zemin))
(if (< zy zemax) (set! zy zemax))
(cons sx zy))
(graph-startpkt fkt (cons (+ sx 1) zy))))))
(cons sx zy))))]
[graph-zeichnen
(lambda (fkt pkt)
(let ([sx (car pkt)] [zy (cdr pkt)])
(if (< sx (+ spbis 1))
(with-failure-continuation
(lambda (error-record error-k)
(graph-zeichnen fkt (graph-startpkt fkt (cons (+ sx 1) zy))))
(lambda ()
(let* ([sx-neu (+ sx 1)]
[x (+ (* (/ (- sx-neu spmin) (- spmax spmin)) (- xo xu)) xu)] [y (fkt x)])
(if (real? y)
(let ([zy-neu (+ (* (/ (- y yu) (- yo yu)) (- zemax zemin)) zemin)])
(if (> zy-neu zemin) (set! zy-neu zemin))
(if (< zy-neu zemax) (set! zy-neu zemax))
(if (and (< zy zemin) (> zy zemax) (< zy-neu zemin) (> zy-neu zemax))
(line canvas sx zy sx-neu zy-neu))
(graph-zeichnen fkt (cons sx-neu zy-neu)))
(graph-zeichnen fkt (graph-startpkt fkt (cons (+ sx 1) zy))))))))))]
[graph-list-zeichnen
(lambda (fkt-ls anz color-ls)
(if (not (null? fkt-ls))
(let ([fkt (eval (car fkt-ls))] [color (n-tes-mod color-ls anz)])
(if (list? color)
(apply setcolor canvas color)
(apply setcolor canvas (eval color)))
(graph-zeichnen fkt (graph-startpkt fkt (cons spvon z0)))
(graph-list-zeichnen (cdr fkt-ls) (+ anz 1) color-ls))))])
(if (list? koordcolor)
(apply setcolor canvas koordcolor)
(apply setcolor canvas (eval koordcolor)))
(if (< s0 spmin) (set! s0 spmin)) (if (> s0 spmax) (set! s0 spmax))
(if (> z0 zemin) (set! z0 zemin)) (if (< z0 zemax) (set! z0 zemax))
(line canvas spmin z0 spend z0)
(line canvas s0 zemin s0 zeend)
(line canvas (- spend 6) (- z0 2) spend z0)
(line canvas (- spend 6) (+ z0 2) spend z0)
(line canvas (- spend 6) (- z0 2) (- spend 6) (+ z0 2))
(line canvas (- s0 2) (+ zeend 6) s0 zeend)
(line canvas (+ s0 2) (+ zeend 6) s0 zeend)
(line canvas (- s0 2) (+ zeend 6) (+ s0 2) (+ zeend 6))
(set! mx (anpassen mx)) (set! teilerx (teilen mx))
(set! stepx (/ (- spmax spmin) mx teilerx))
(set! my (anpassen my)) (set! teilery (teilen my))
(set! stepy (/ (- zemin zemax) my teilery))
(x-teilen-positiv (+ s00 stepx)) (x-teilen-negativ (- s00 stepx))
(y-teilen-positiv (- z00 stepy)) (y-teilen-negativ (+ z00 stepy))
(set! ex (/ ex mx teilerx))
(set! faktorx (faktor
(* 3/2 textgroesse (string-length (number->string (genaurunden ex genau))))
stepx))
(set! ey (/ ey my teilery))
(set! faktory (faktor
(* 3/2 textgroesse)
stepy))
(if (<= (abs yu) (abs yo))
(begin
(x-bezeichnen-positiv (+ s00 (* stepx faktorx)) (* ex faktorx) #t)
(x-bezeichnen-negativ (- s00 (* stepx faktorx)) (- (* ex faktorx)) #t))
(begin
(x-bezeichnen-positiv (+ s00 (* stepx faktorx)) (* ex faktorx) #f)
(x-bezeichnen-negativ (- s00 (* stepx faktorx)) (- (* ex faktorx)) #f)))
(if (<= (abs xu) (abs xo))
(begin
(y-bezeichnen-positiv (- z00 (* stepy faktory)) (* ey faktory) #t)
(y-bezeichnen-negativ (+ z00 (* stepy faktory)) (- (* ey faktory)) #t))
(begin
(y-bezeichnen-positiv (- z00 (* stepy faktory)) (* ey faktory) #f)
(y-bezeichnen-negativ (+ z00 (* stepy faktory)) (- (* ey faktory)) #f)))
(x-achse-bezeichnen bez)
(y-achse-bezeichnen bez)
(if (< spvon spmin) (set! spvon spmin)) (if (> spbis spmax) (set! spbis spmax))
(graph-list-zeichnen fkt-list 0 color-list)
'ok))))
</code>
<code id="53ea65952df22">
(define runden
(lambda (x eps)
(let ([anzahl (inexact->exact (round (- (/ (log eps) (log 10)))))])
(exact->inexact (/ (round (* x (expt 10 anzahl))) (expt 10 anzahl))))))
(define termeingabe
(lambda ()
(let ([proz (read (open-input-string (inputDialog "Prozedur (Lambda-Term):")))])
(with-failure-continuation
(lambda (error-record error-k) #f)
(lambda () (eval proz))))))
(define fkt-list
(lambda (fkt x)
(letrec ([suchen
(lambda (x anz ls)
(cond
[(< anz 1) (reverse ls)]
[(< (abs (fkt x)) (expt 10 (- genau))) (reverse ls)]
[else
(let* ([m ((ableitung fkt) x)] [y (fkt x)] [n (- y (* m x))])
(suchen (- x (/ y m)) (- anz 1) (cons `(lambda (x) (+ (* ,m x) ,n)) ls)))]))])
(suchen x maxanz (list fkt)))))
(define nullstellenausgabe
(lambda (fkt x)
(let ([x0 (nst fkt x)])
(display "Nullstelle: ")
(if x0
(display (number->string (runden x0 eps)))
(display "--> keine"))
(newline)
'ok)))
(define nullstelle
(lambda (xu xo yu yo x-start)
(let ([fkt (termeingabe)])
(if fkt
(let* ([ls (fkt-list fkt x-start)] [anz (length ls)])
(schaubild ls xu xo yu yo xu xo black (rgb-list anz))
(nullstellenausgabe fkt x-start))
(begin
(display "--> Eingabefehler!") (newline)
'ok)))))
</code>
'''Hinweise:'''
* Die Prozedur <tt>rgb-list</tt> liefert eine beliebige Anzahl von aufeinander abgestimmten Farbwerten.
* Die zu untersuchende mathematische Funktion wird als Lambda-Term in einem ''Interaktionsfenster'' eingegeben, z.B.:
(lambda (x) (+ (expt x 3) (* -2 (expt x 2)) (* -11 x) 12))
<canvas w="610" h="460"></canvas>
<run id="53ea65952e175">
(nullstelle -10 10 -25 25 3)
</run>
Zurück zum [[Grundlagen der funktionsorientierten Programmierung mit SCHEME/Der Algorithmusbegriff|Algorithmusbegriff]].