Bei einer Menge von Punkten in möchte ich berechnen genau. ist das Lagrange-Polynom in Bezug auf die Punkte mit als Knoten, dh Da dies ein Polynom vom Grad , könnte ich jede alte Gaußsche Quadratur von ausreichendem Grad verwenden. Dies funktioniert gut, wenn nicht zu groß ist, führt jedoch zu Ergebnissen, die durch Rundungsfehler für großes fehlerhaft sind . [ - 1 , 1 ] ∫ 1 - 1 L i ( x )L i x j x i L i ( x ) = ∏ j ≠ i x - x j
nnn
Irgendeine Idee, wie man diese vermeidet?
quadrature
integration
Nico Schlömer
quelle
quelle
Antworten:
Die Berechnung von für die Lagrange-Polynome die auf einem beliebigen Gitter kann von der durchgeführt werden folgende zwei Schritte:∫1−1Lk(x)dx Lk xk,k=0,…,n
Berechnen Sie die Clenshaw-Curtis-Quadraturgewichte auf dem Chebyshev-Extrema-Gitter für : wobei für , andernfalls und für oder , andernfalls . Siehe das Papier von Waldvogel (2006) für weitere Details.wcck yk k=0,…,n yk=cos(kπn)wcck=ckn(1−∑j=1⌊n/2⌋bj4j2−1cos(2πjkn)) bk:=1 k=n/2 bk:=2 ck:=1 k=0 k=n ck:=2
Transformiere die Gewichte über die Transformationsmatrix in das beliebige Gitter , um die gesuchten Gewichte , wobeiwcck xk,k=0,…,n M wk wk=∑jMkjwccj Mjk = Lj(yk).
Im Prinzip ist dies nur eine Clenshaw-Curtis-Quadratur mit Funktionswerten auf dem beliebigen Gitter , die jedoch durch Basistransformation erhalten wird (für eine allgemeine Referenz zu Clenshaw-Curtis siehe z. B. das Trefethen-Papier).xk
Der Algorithmus scheint ziemlich stabil zu sein, insbesondere im Vergleich zum Vandermonde-Ansatz, wie er in der Antwort von @Kirill angegeben ist : Obwohl er denselben Vorstellungen folgt - generieren Sie die Quadraturgewichte auf bekannter Basis und transformieren Sie sie dann in das neue Gitter - dies hätte erwartet werden können, da die Transformation in Bezug auf die Vandermonde-Matrix normalerweise stark schlecht konditioniert ist.
Beispiel: Erzeugung von Legendre-Lobatto-Quadraturgewichten
Wir betrachten das Beispiel der Legendere-Lobatto-Quadraturregel und vergleichen die Genauigkeit mit dem monomialen Ansatz. Als Referenz verwenden wir die Quadraturgewichte die vom Golub-Welsch-Algorithmus für verschiedene und berechnen den kumulierten Fehler Hier ist das Ergebnis: Man stellt fest, dass die Clenshaw-Curtis-Quadraturgewichte über den betrachteten Bereich von Gitterpunkten perfekt stabil sind und die Legendre-Gewichte bis zur Maschinengenauigkeit ( reproduzieren ).wLegk n ϵn=∑k=1n(wk−wLegk)2 ϵ√∼1015
Beispiel: Erzeugung von Newton-Cotes-Quadraturformeln
Wir betrachten die Erzeugung der Newton-Cotes-Quadraturformel auf gleich beabstandeten Gittern. Wiederum erwartet man eine schlechte Konditionierung, da kurz gesagt für die Polynominterpolation Gitter mit gleichem Abstand baaad sind.
Im folgenden Bild habe ich die absolute Summe der Gewichte berechnet .∑i|wi|/N
Bis zu beispielsweise 50 Gitterpunkte stimmen die Ergebnisse des Monomial- und des Clenshaw-Curtis-Ansatzes überein. Danach wird Clenshaw-Curtis besser - für das, was es wert ist. Eine direkte Interpretation ist, dass das gleichmäßig verteilte Gitter alles für beispielsweise ruiniert . Bei etwa schlägt der Zustand der Vandermonde-Matrix jedoch zurück und führt zu einem noch schlechteren Ergebnis.n>10 n=50
Beispiel: Guass-Patterson-Quadratur
Dieses Beispiel ist @ NicoSchlömer zu verdanken. Ich kannte diese Regeln bisher nicht, also habe ich die Abszissen aus dieser Implementierung genommen und sowohl den Vandermonde- als auch den transformierten Clenshaw-Curtis-Ansatz angewendet (wobei der Vandermonde-Ansatz wie oben den Björk-Pereyra-Algorithmus verwendet).
Wie im Kommentar vorgeschlagen, berechnete ich dann den Fehler der Integration einer konstanten Funktion durch mit dem folgendes Ergebnis:ϵ=1n∣∣∣2−∑i=1nwi∣∣∣,
Aus diesem Bild geht hervor, dass der transformierte Clenshaw-Curtis-Ansatz weitaus effizienter ist als der Vandermonde-Ansatz (zumindest in der Arithmetik mit endlicher Genauigkeit). Trotzdem bricht Clenshaw-Curtis ab Index 7 zusammen, daher sollten andere Methoden angewendet werden.
quelle
Sie können dies mit dem Björck-Pereyra-Algorithmus zur Lösung von Vandermonde-Systemen bewerten, da Sie bewertenb⊤V−1 b=(2,0,23,0,25,0,…)
Ausgabe:
Wenn dies
X
nicht wie in diesem Test positiv ist, scheinen die relativen Fehler in der gleichen Größenordnung zu liegen wie bei einer regulären linearen Lösung.quelle
V.main(32)
auf meinem Laptop in etwa einer Sekunde eine sinnvolle Antwort erzeugt (mit nur wenig Speicher). Die Zahlen sind nicht einmal so groß, der größte Zähler hat 54 Ziffern, also vermute ich, dass etwas anderes für Sie schief geht. Kannst du einen Kern posten, weil ich gespannt bin, wie es fehlschlägt?Float64
fürd
: check with@show typeof(d)
. Lassen Sie mich wissen, wenn Sie weitere Probleme damit finden.Berechnen Sie zuerst die Produkte der Nominatoren und Nenner und teilen Sie sie dann einmal. Die beiden Produkte sollten in der gleichen Größenordnung liegen, damit keine signifikanten Rundungsfehler auftreten. Durch die geringere Anzahl von Gleitkommaberechnungen erhalten Sie außerdem den zusätzlichen Vorteil einer höheren Geschwindigkeit.
quelle