Numerische Lösungen mathematischer Probleme

Wir stellen Ihnen eine kleine Sammlung von mathematischen Unterprogrammen vor, die in eigenen Programmen benötigt werden, aber schwer zu programmieren sind.
In dem ersten Programm (Listing 1) wird näherungsweise die erste Ableitung einer Funktion an einer bestimmten Stelle mit Hilfe des Differenzenquotienten bestimmt. Der Differenzenquotient ist wie folgt definiert:

\( f'(z) = \frac{f(x) - f(x_0)}{\,x - x_0\,} \quad\text{mit}\quad x_0 \le z \le x\)

Dieser Differenzenquotient (im Programm Zeile 200) wird etwas abgewandelt, so daß man den Abstand zwischen x und x0 immer kleiner werden läßt Dies geschieht, indem man statt x den Wert xo+h nimmt. Der Wert h wird dann immer weiter halbiert (Zeile 230), bis eine vorgegebene Genauigkeit erreicht ist (Zeile 220). Als erste Ableitung der bestimmten Stelle ist dann der letzte ausgedruckte Wert zu nehmen. Die erste Ableitung gibt übrigens die Steigung der Tangenten an, die an der bestimmten Stelle an den Graphen der Funktion angelegt wird.

Diese Tangente, beziehungsweise auch weitere Tangenten werden bei dem zweiten und dritten Programm (Listing 2 und 3) benötigt. Bei dem zweiten Programm wird mit Hilfe von Tangenten die Nullstelle(n) (die Stelle, wo der Graph der Funktion die x-Achse schneidet) einer Funktion bestimmt.

Die Iterationsvorschrift (Zeile 160) lautet:

\( x_{i+1} = x_i - \frac{y_i}{y_i'} = x_i - \frac{f(x_i)}{f'(x_i)}\)

Dieses Verfahren stammt von Newton. Im dritten Programm (Listing 3) ist die Vorschrift ähnlich, nur braucht man die erste Ableitung nicht wissen, was oft von Vorteil ist. Die erste Ableitung wird hier in etwa so bestimmt, wie es in Programm 1 beschrieben ist. Da nicht nur Ableitungen und Nullstellen von Interesse sind, sondern auch die oben erwähnten Integrale, behandeln Programm 4 und 5 (Listing 4 und 5) diese Probleme.

Im vierten Programm (Listing 4) ist die Kepler’sche Faßregel umgesetzt. Möchte man das Integral einer Funktion (Zeile 240) wissen, so kann man es mit der Kepler’schen Faßregel annähern. Die Faßregel lautet:

\( \int_a^b f(x)\,dx \;\approx\; (b-a)\,\Bigl(\tfrac{f(a)}{6} \;+\; \tfrac{2}{3}\,f\!\bigl(\tfrac{a+b}{2}\bigr) \;+\; \tfrac{f(b)}{6}\Bigr)\)

Mit dieser noch recht einfachen Formel läßt sich das Integral schnell näherungsweise bestimmen.

Die Formel im Programm 5 (Listing 5) stammt von Simpson, sie ist ein wenig komplizierter:

\( \int_a^b f(x)\,dx \;\approx\; \frac{h}{3}\,\Bigl[f(a) + f(b) + 4\sum_{k=0}^{n-1} f\!\bigl(a + (2k+1)h\bigr) + 2\sum_{k=1}^{n-1} f\!\bigl(a + 2kh\bigr)\Bigr]\)
mit \( h = \frac{b-a}{2n}\)

Die beiden Summen sind in den Zeilen210bis 240 und 250 bis 280 umgesetzt.

Ein Problem aus einem etwas anderen Bereich bezieht sich auf gewöhnliche Differentialgleichungen. Bei dem Eulerschen Polygonzugverfahren in Programm 6 (Listing 6) wird der Graph, der die Lösung der gewöhnlichen Differentialgleichung beschreibt, in Form von kleinen Geradenstückchen genähert. Diese kleinen Geradenstückchen bilden dann den Polygonzug. Die einzelnen Punkte zu dieser Näherung werden auf einem vorgegebenen Bereich bezüglich eines bestimmten Anfangswertes vorgegeben.

Das Verfahren von Runge-Kutta in Programm 7 (Listing 7) ist weitaus genauer, aber dafür ist die Formel auch komplizierter (Zeile 200, Zeilen 250 bis 280):
Anfangswert \( y_0 = y(x_0)\)
Schrittfunktion (Schrittweite h):
\( y_{k+1} = y_k + \frac{h}{6}\,(L_{k,1} + 2L_{k,2} + 2L_{k,3} + L_{k,4})\)
\(\text{mit } L_{k,1} = f(x_k, y_k) \;\; \{\text{mit } y' = f(x,y)\}\)
\( L_{k,2} = f\!\bigl(x_k + \tfrac{h}{2},\, y_k + \tfrac{h}{2}L_{k,1}\bigr)\)
\( L_{k,3} = f\!\bigl(x_k + \tfrac{h}{2},\, y_k + \tfrac{h}{2}L_{k,2}\bigr)\)
\( L_{k,4} = f\!\bigl(x_k + h,\, y_k + h\,L_{k,3}\bigr)\)

100 rem differenzenquotient
110 print"differenzenquotient"
120 input"differentiationsstelle";x0
130 input"startstelle           ";x1
140 ifx0=x1then120
150 rem berechnung und ausgabe
160 print" x"," f'"
170 h=x1-x0
180 x=x0+h:gosub270:f0=f
190 x=x0:gosub270:f1=f
200 df=(f0-f1)/h
210 printint((x0+h)*1e5+.5)/1e5,int(df*1e5+.5)/1e5
220 ifh<>(x1-x0)andabs(df-dg)<1e-6then250
230 dg=df:h=h/2
240 goto180
250 print"f'(";x0;")=";df
260 end
270 rem funktion
280 f=x*x*x-2*x*x-2
290 return

100 rem nullstellensuche nach newton
110 print"iteration nach newton"
120 print" x"," f"
130 readx0:rem startwert
140 x=x0:gosub220:f0=f:rem funktionswert
150 x=x0:gosub250:f1=fs:rem 1. ableitung
160 x1=x0-f0/f1
170 printx1,f0
180 ifabs(x1-x0)<1e-6thenend
190 x0=x1
200 goto140
210 rem funktion
220 f=x^2-4
230 return
240 rem 1. ableitung
250 fs=2*x
260 return
270 rem startwert
280 data3

100 rem nullstellensuche - regula-falsi
110 print"regula falsi"
120 readx0,x1:rem startwerte
130 print" x"," f"
140 x=x1:gosub250:f1=f
150 x=x0:gosub250:f0=f
155 printx,f0
160 z=x1-(x1-x0)/(f1-f0)*f1
180 x0=x1:x1=z
190 ifabs(x1-x0)<1e-8then210
200 goto140
210 rem ausgabe
220 print"nullstelle=";z
230 end
240 rem funktion
250 f=x^2-4
260 return
270 rem startwerte
280 data3,4

100 rem integration nach kepler
110 print"keplersche fassregel"
120 input"interervallgrenzen";a,b
130 ifb<=athenprint"falsche eingabe!":goto120
140 rem berechnung des genaeherten integralwertes
150 x=a:gosub240
160 i=i+f/6
170 x=b:gosub240
180 i=i+f/6
190 x=(a+b)/2:gosub240
200 i=i+f*2/3:i=i*(b-a)
210 rem ausgabe
220 print"ergebnis=";i
230 end
240 f=x*x*2
250 return

100 rem integration nach simpson
110 print"simpsonsche regel"
120 input"interervallgrenzen";a,b
130 ifb<=athenprint"falsche eingabe!":goto120
140 input"teilintervalle";n
150 h=(b-a)/(2*n):rem schrittweite
160 rem berechnung des genaeherten integralwertes
170 x=a:gosub330
180 i=i+f
190 x=b:gosub330
200 i=i+f
210 forj=0ton-1
220 x=a+(2*j+1)*h:gosub330
230 i=i+4*f
240 nextj
250 forj=1ton-1
260 x=a+2*j*h:gosub330
270 i=i+2*f
280 nextj
290 i=i*h/3
300 rem ausgabe
310 print"ergebnis=";i
320 end
330 f=x*x*2
340 return

100 rem eulersches polygonverfahren
110 :
130 print"{clr}{rvon}     eulersches polygonzugverfahren     {rvof}"
140 print"dieses verfahren hat keine fehler-"
150 print"abschaetzung. numerisch ist dieses "
160 print"verfahren recht schnell, aber bei zu"
170 print"grosser schrittweite sehr ungenau, wobei"
180 print"{up}der fehler mit wachsendem x immer "
190 print"groesser wird. schrittweiten von 10^-2"
200 print"und kleiner sind schon einigermassen"
210 print"genau.{down}""
220 :
230 input"{down}anfangswerte : x,y ";x1,y1
240 input"{down}intervallobergrenze";b:a=x1
250 input"{down}schrittweite       ";h
260 :
280 print"{clr}y'=y*x-2*x{down}"
290 deffnf(x)=int(x*10e6+.5)/10e6
300 :
310 y=y1
320 print" x   ","  y"
330 forx=atobsteph
340 gosub410:m=f:rem steigung der naeherungsgeraden
350 t=y-m*(x+a):rem verschub der geraden
360 printfnf(x),fnf(y)
370 y=m*(a+x+h)+t:rem naechster y-wert
380 nextx
390 end
400 :
410 f=y*x-2*x:rem differentialgleichung
420 return

100 print"{clr}{rvon} runge-kutta-verfahren ohne lok. fehler {rvof}"
110 print"dgl am programmende eingeben!"
120 input"{down}schrittweite  ";h
130 input"{down}anfangswert   ";y0
140 input"{down}interv.grenzen";a,b
150 deffnr(x)=int(x*1e4+.5)/1e4
160 print"{down}{down} x"," y"
170 yk=y0
180 xk=a
190 gosub250
200 yl=yk+(h/6)*(l1+2*l2+2*l3+l4)
210 printfnr(xk),yk
220 xk=xk+h:yk=yl
230 ifxk<=b+hthen190
240 end
250 x=xk:y=yk:gosub300:l1=f
260 x=xk+h/2:y=yk+(h/2)*l1:gosub300:l2=f
270 x=xk+h/2:y=yk+(h/2)*l2:gosub300:l3=f
280 x=xk+h:y=yk+h*l3:gosub300:l4=f
290 return
300 f=y*x-2*x  :rem <=> y'=... _dgl
310 return