Quadraturalgorithmen

26.10.2020, 11:25

Finn_

Quadraturalgorithmen

Diese Woche wollte ich mich mal um einen besseren Quadratur-Algorithmus kümmern, wofür irgendwie eine große Notwendigkeit besteht. Bisher habe ich für meinen Funktionenplotter und andere Software die Gauss-Legendre-Quadratur mit 64 Stützstellen genutzt. Die Berechnung geht so schnell, dass das Plotten von Integralfunktionen kein Problem darstellt.

Nun verhält es sich ja allerdings so, dass diese bei unstetigen Funktionen bzw. Funktionen mit hoher Krümmung nicht ausreichend genau ist. Natürlich gibt es einen Parameter n mit Vorgabewert n=1, der das Intervall in n Teile zerlegt. Die Berechnung wird dann meistens ausreichend genau, allerdings muss man dafür ein wenig herumspielen und der Berechnungsaufwand steigt entsprechend.

Gesucht ist schlicht, was man unter eierlegende Wollmilchsau versteht. Das will heißen, es sollte eine adaptive Berechnung erfolgen, bei der man n in den meisten Fällen nicht mehr anfassen braucht.

Ich möchte kein perfektes Verfahren (dafür scheint ohnehin Intervall-Arithmetik notwendig zu sein, das wird kompliziert), sondern etwas das mir im Alltag quick and dirty in den allermeisten Fällen das richtige Resultat liefert. Außerdem sollte die Implementation möglichst kompakt sein, damit man sie unter widrigen Umständen schnell auf andere Systeme und Programmiersprachen portieren kann.

Meine erste Überlegung ist folgender Algorithmus.

code:

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:

GL32 = [
[-0.9972638618494816, 0.007018610009470141],
[-0.9856115115452684, 0.01627439473090563],
[-0.9647622555875064, 0.02539206530926208],
[-0.9349060759377397, 0.03427386291302148],
[-0.8963211557660522, 0.04283589802222668],
[-0.8493676137325699, 0.05099805926237622],
[-0.7944837959679425, 0.05868409347853551],
[-0.7321821187402897, 0.06582222277636184],
[-0.6630442669302152, 0.07234579410884850],
[-0.5877157572407623, 0.07819389578707028],
[-0.5068999089322295, 0.08331192422694673],
[-0.4213512761306354, 0.08765209300440380],
[-0.3318686022821277, 0.09117387869576390],
[-0.2392873622521371, 0.09384439908080457],
[-0.1444719615827965, 0.09563872007927489],
[-0.04830766568773831,0.09654008851472778],
[ 0.04830766568773831,0.09654008851472778],
[ 0.1444719615827965, 0.09563872007927489],
[ 0.2392873622521371, 0.09384439908080457],
[ 0.3318686022821277, 0.09117387869576390],
[ 0.4213512761306354, 0.08765209300440380],
[ 0.5068999089322295, 0.08331192422694673],
[ 0.5877157572407623, 0.07819389578707028],
[ 0.6630442669302152, 0.07234579410884850],
[ 0.7321821187402897, 0.06582222277636184],
[ 0.7944837959679425, 0.05868409347853551],
[ 0.8493676137325699, 0.05099805926237622],
[ 0.8963211557660522, 0.04283589802222668],
[ 0.9349060759377397, 0.03427386291302148],
[ 0.9647622555875064, 0.02539206530926208],
[ 0.9856115115452684, 0.01627439473090563],
[ 0.9972638618494816, 0.007018610009470141]
]

# Gauss-Legendre-Quadratur mit 32 Stützstellen
def gauss(f,a,b):
    p = 0.5*(b-a); q = p+a
    return p*sum(w*f(p*x+q) for x,w in GL32)

# Rekursive Verfeinerung bis Krümmung niedrig genug ist
def quad_rec(a,b,f,depth):
    if depth==0: return gauss(f,a,b)
    h = 0.1*(b-a)
    y = [f(a+k*h) for k in range(11)]
    dmax = max(abs(0.5*(y[k+2]+y[k])-y[k+1]) for k in range(9))
    if dmax < 0.001:
        return gauss(f,a,b)
    else:
        return sum(quad_rec(a+k*h,a+(k+1)*h,f,depth-1)
            for k in range(10))

def quad(a,b,f,n=1):
   h = (b-a)/n
   return sum(quad_rec(a+k*h,a+(k+1)*h,f,10) for k in range(n))

Betrachten wir dieses Testintegral aus den Folien "Fast and rigorous numerical integration" von Fredrik Johansson:

$\begin{eqnarray*} \int_0^8 (e^x-\lfloor e^x\rfloor)\sin(x+e^x)\,\mathrm dx. \end{eqnarray*}$

Der korrekte Wert ist
0.098651704478365206119658...

[attach]52047[/attach]

code:

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:

from math import floor, exp, sin

count = 0
def f(x):
    global count
    count += 1
    return (exp(x)-floor(exp(x)))*sin(x+exp(x))

Delta = abs(quad(0,8,f)-0.098651704478365206119658)

print("Abweichung: {:.2E}".format(Delta))
print("Anzahl der Funktionsaufrufe:",count)

# Ausgabe:
# Abweichung: 9.59E-11
# Anzahl der Funktionsaufrufe: 6980471

Tja, die Genauigkeit genügt mir, aber zum Zeichnen der Integralfunktion ist das zu langsam.

26.10.2020, 12:22

HAL 9000

Auf diesen Beitrag antworten »

Zitat:

Original von Finn_
# Anzahl der Funktionsaufrufe: 6980471

Ich schätze mal, es dürfte schwierig sein, DIE eierlegende Wollmichsau in diesem Gebíet zu finden, wenn sie noch dazu solch exotisches Futter wie den obigen (zumindest zum Ende zu) wild oszillierenden Integranden $\begin{eqnarray*} \left(e^x-\lfloor e^x\rfloor\right)\sin\left(x+e^x\right) \end{eqnarray*}$ schnell (!) verdauen können muss.

Die ABSOLUTE Schranke 0.001 in der Abbruchbedingung ist mir nicht geheuer: D.h., wendet man den Algorithmus statt auf $\begin{eqnarray*} f \end{eqnarray*}$ z.B. auf $\begin{eqnarray*} 1000\cdot f \end{eqnarray*}$ an, wird länger (und mutmaßlich genauer) gerechnet? Eine solche Entscheidung sollte doch eher skalenunabhängig gefällt werden - ist zumindest meine Meinung.

Warum die Verfeinerung überhaupt in Zehnerpotenzschritten geschieht, wäre auch noch zu hinterfragen. Ich würde Zweierpotenzschritte als natürlicher empfinden, wobei man am besten gleich noch die Berechnung der Funktionswerte einspart, die man im vorigen Durchlauf schon hatte - ist besonders bei "zeitteuren" Funktionswertberechnungen anratsam.

Neue Frage »

Antworten »

Quadraturalgorithmen

Verwandte Themen