Berechnung einer Quadraturformel zur Integration von multivariaten Polynomen über sphärische Dreiecke

Judith Beckmann
judith.beckmann@freenet.de

Beschreibung des Problems

Die Einheitssphäre im euklidischen Raum (𝕊^2)={x∈(ℝ^3): ||x(Subscript[||, 2])=1} kann 
Modell der Erdoberfläche in der Meteorologie und Geophysik sein oder als Oberfläche von Organen in der medizinischen Modellierung dienen. Satelliten, Flugzeuge, GPS u.a. liefern Daten über die Struktur oder physikalische Phänomene an verschiedenen Punkten der Erdoberfläche. In der Medizin erhalten wir Daten durch bildgebende Verfahren wie z.B. Ultraschall und MRT. Es liegen also in vielen praktischen Anwendungen große Datenmengen nur in Form von Funktionswerten vor, deren Lokalisation oft nicht beeinflusst werden kann, sog. Scattered Data.

In vielen Anwendungen geht es darum, eine glatte Funktion auf der Sphäre zu finden, die die gegebene Datenmenge interpoliert oder approximiert, oder darum, Integral- und 
Differentialgleichungen zu lösen. Dazu brauchen wir Quadraturformeln zur Integration auf der Einheitssphäre. In [8] wird gezeigt, dass Quadraturformeln mit positiven Gewichten, die exakt für Polynome sind, auf der Sphäre existieren. Die Existenz von ähnlichen Quadraturformeln auf sphärischen Kappen oder kompakten Untermengen der Einheitssphäre wird in [9] und [10] bewiesen. 

Die Sphäre und jede kompakte Untermenge auf der Einheitssphäre kann in sphärische Dreiecke aufgeteilt, d.h. trianguliert werden. Bei großen Datenmengen können so Speicherplatzprobleme gelöst werden.  
Ist die den Daten zugrundeliegende Funktion nur in einem kleinen Gebiet auf der Sphäre gut lokalisiert, dann ist es für die Approximation des Integrals durch eine Quadraturformel sicherlich von Vorteil, in dem entsprechenden Gebiet möglichst viele Funktionswerte zu kennen, während wir in den Restgebieten weniger Informationen über die Funktion benötigen. Wenn die Quadraturgewichte speziell für das Gebiet, in dem die Funktion gut lokalisiert ist, berechnet werden, kann das Konvergenzverhalten der Quadratur möglicherweise schon bei kleinen Polynomgraden günstig sein.

In [2] wird die Konstruktion von Quadraturformeln konkret über sphärische Dreiecke beschrieben nach einer Idee von Le Gia und Mhaskar in [7]. Da wir bisher über keine zu den Spherical Harmonics vergleichbaren orthogonale Polynome über kompakte Untermengen verfügen, müssen wir Integrale auf diesen Untermengen für Spherical Harmonics exakt berechnen können.


Sei (Subscript[ℍ, n]^2) der Raum der Spherical Harmonics vom Höchstgrad n, dim((Subscript[ℍ, n]^2))=(((n+1))^2). Eine Orthonormalbasis des (Subscript[ℍ, n]^2) bezüglich des (L^2)((𝕊^2))-Skalarproduktes sei 
                   (Subscript[Y, q,m])(θ,ρ):=(Sqrt[((2q+1)/(4π))])(Subscript[P, q]^(|m|))(cos θ) (e^(imρ)), q=0,1,...,n, m=-q,...,q. 
Es seien die Datenpaare ((Subscript[y, k]),f((Subscript[y, k]))) auf einem sphärischen Dreieck Δ((Subscript[V, 1]),(Subscript[V, 2]), (Subscript[V, 3])) mit den Ecken (Subscript[V, 1]),(Subscript[V, 2]),(Subscript[V, 3 ])∈  (𝕊^2), (Subscript[y, k]) ∈  Δ((Subscript[V, 1]),(Subscript[V, 2]),(Subscript[V, 3])) ⊂ (𝕊^2), f((Subscript[y, k])) ∈ ℂ, k=1,2,...,M, M > (((n+1))^(2 ))gegeben.
Gesucht sind Quadraturgewichte (Subscript[λ, k]) zu den Knoten (Subscript[y, k]), so dass
                                        (∫_(Δ((Subscript[V, 1]), (Subscript[V, 2]), (Subscript[V, 3]))))P(x)  dξ(x)=(∑_(k=1)^M)(Subscript[λ, k])P((Subscript[y, k])) für alle P∈ (Subscript[ℍ, n]^2).

Zu lösen ist ein unterbestimmtes Lineares Gleichungssystem Yα=b mit 
(Subscript[(Y), l,k])=(Subscript[Y, l])((Subscript[y, k])), l=1,2,...,(((n+1))^2), k=1,2,...,M, und (Subscript[(b), l]):=(∫_(Δ((Subscript[V, 1]), (Subscript[V, 2]), (Subscript[V, 3]))))(Subscript[Y, l])(x)  dξ(x).

Hat Y Rang (((n+1))^2), was wir für M ≫ (((n+1))^2) voraussetzen können, dann erhalten wir als eindeutige Lösung des Kleinste-Quadrate-Problems
                                                               (min_(α∈(ℂ^M)))||(W^(-1/2))α(Subscript[||, 2]^2), Yα=b, 
W ist eine Diagonalmatrix mit positiven reellen Einträgen auf der Diagonalen, mit der Normalengleichung zweiter Art die Gewichte λ=((Subscript[λ, k])), k=1,2,...,M. 

 Dabei bleiben die Integrale (Subscript[(b), l])=(∫_(Δ((Subscript[V, 1]), (Subscript[V, 2]), (Subscript[V, 3]))))(Subscript[Y, l])(x)  dξ(x) zu lösen.
 
 
 
 
In [1] konstruieren wir eine Quadraturformel, mit der wir Integrale                                                                     (∫_(Δ((Subscript[V, 1]), (Subscript[V, 2]), (Subscript[V, 3]))))(Subscript[Y, n])(x)  dξ(x), (Subscript[Y, n])∈ (Subscript[ℍ, n]^2),
exakt berechnen können. 
Mit unserem Algorithmus bestimmen wir Gewichte (Subscript[w, r,s]) und Knoten ((Subscript[χ, r]),(Subscript[ψ, s])), 1⩽|r|⩽n+2, 1⩽|s|⩽n+1, so dass
           (∫_(Δ((Subscript[V, 1]), (Subscript[V, 2]), (Subscript[V, 3]))))(Subscript[Y, n])(x)  dξ(x)=(∑_(1⩽|r|⩽n+2))(∑_(1⩽|s|⩽n+1))(Subscript[w, r,s])(Subscript[Y, n])((Subscript[χ, r]),(Subscript[ψ, s])) für alle (Subscript[Y, n])∈(Subscript[ℍ, n]^2).
Mit der Implementierung der Funktionen im Mathematica Package GaussianQuadratureSphericalTriangle.m können wir die einzelnen Berechnungsschritte numerisch realisieren.

Anmerkungen zum Algorithmus und zur Implementierung:

Ein beliebiges sphärisches Dreieck Δ((Subscript[V, 1]),(Subscript[V, 2]),(Subscript[V, 3])) wird, wenn nötig, an den Nordpol verschoben und gegebenenfalls in drei gleichschenklige Dreiecke (Subscript[Δ, j]), j=1,2,3, zerlegt (Funktion Translation), z.B. wenn es selber nicht gleichschenklig ist.

Berechnet werden die Quadraturgewichte und -knoten für jedes gleichschenklige Dreieck (Subscript[Δ, j]) , wobei die Ecke gegenüber der Basis der Nordpol (Subscript[n, 0])=(0,0,1) ist. Für den Winkel γ gegenüber der Basis (vertex angle) gilt 0<γ<π. Das Dreieck liegt in einer sphärischen Kappe (Subscript[𝕊, ρ]^2):={x∈(𝕊^2):<x,(Subscript[n, 0])(Subscript[>, 2]) ⩾ cos ρ} mit Radius ρ, 0<ρ<π/2. Der Spezialfall ρ=π/2 wird in der Implementierung des Packages berücksichtigt. Aus den drei Quadraturen ergibt sich die Quadraturformel für das Dreieck Δ((Subscript[V, 1]),(Subscript[V, 2]),(Subscript[V, 3])) (Funktion GaussQuadratureST mit Hilfsfunktionen).

Die Funktionen CIRCUM, LEFT und SCOORD sind Übertragungen aus Renkas Fortran 77 Software-Paket STRIPACK (siehe [11]). Sie bestimmen den Mittelpunkt 
des Dreiecksumkreises, die Position eines Punktes auf der Sphäre bzgl. eines Kreisbogens bzw. liefern die sphärischen Koordinaten für einen Punkt in kartesischen Koordinaten.

In [1] werden im 5. Kapitel die einzelnen Berechnungsschritte zusammenfassend dargestellt: 
Wir benötigen Quadraturformeln über Kreisbögen bzgl. verschiedener Gewichtsfunktionen und Kerne für beliebige trigonometrische Polynome über Kreisbögen. Die Quadraturen über Kreisbögen beruhen auf den bekannten Gauss-Quadraturen.
Zur Berechnung der benötigten Gauss-Quadraturgewichte und -knoten und Werte der Kerne nutzen wir die Ideen von Gautschi in [3-6] zur Diskretisierung von Gewichtsfunktionen. Die Berechnungsschritte nach den Ideen von Gautschi bezeichnen wir als Gautschi-Algorithmus. 

Der Gautschi-Algorithmus liefert approximierte Dreiterm-Rekursions-Koeffizienten. Dazu benötigt er Gauss-Jacobi-Quadraturen zu Gewichtsfunktionen (1-x)^r, r=-1/2 oder r=1/2 auf dem Intervall (-1,1); die Knoten sind Nullstellen der entsprechenden orthogonalen Jacobi-Polynome und liegen bis zum Polynomgrad 981 (r=-1/2) bzw. 648 (r=1/2) im Verzeichnis Zeros vorberechnet vor (bis Polynomgrad 300 bzw. 100 im Verzeichnis Zeros_part).

Der Gautschi-Algorithmus liegt in Mathematica realisiert in der Funktion mcdismod vor und verwendet wahlweise einen diskretisierten Stieltjes-Algorithmus oder einen modifizierten Lanczos-Algorithmus. Für das Abbruchkriterium benötigen wir, jeweils für die Dreiecke (Subscript[Δ, 1]), (Subscript[Δ, 2]),  (Subscript[Δ, 3]), drei Fehlerschranken (Subscript[ϵ, 1]), (Subscript[ϵ, 2]), (Subscript[ϵ, 3 ]) (wird Δ((Subscript[V, 1]),(Subscript[V, 2]),(Subscript[V, 3])) nicht zerlegt, dann ist nur (Subscript[ϵ, 1])relevant). Mit diesen lässt sich der Approximationsfehler der Dreiterm-Rekursions-Koeffizienten und damit die Genauigkeit der Quadraturformel beeinflussen.
Die Funktionen Jacobi, lanczos, stieltjes sind Übertragungen aus Gautschis Matlab Software-Paket OPQ (siehe [5] und [6]), ebenso die Funktion mcdismod mit Anpassungen an unser Problem.

Nähert sich in einem der Teildreiecke (Subscript[Δ, j]) der Winkel (Subscript[γ, j]) dem Winkelmaß π an, dann steigt die Berechnungszeit deutlich an. In unserer Implementierung werden die Quadraturgewichte für ein Dreieck (Subscript[Δ, j]) nicht berechnet, falls |π-(Subscript[γ, j])| ⩽ 10*eps0value. Der Wert von eps0value kann gewählt werden, ist ansonsten gleich (10^(-16)).

Die Verkleinerung des Radius ρ und Zunahme des Polynomgrades n erhöhen die Berechnungszeit. Ein deutlich verkleinerter Radius führt zu einem spürbaren Präzisionsverlust. Mit Zunahme des Polynomgrades n nimmt die Präzision linear ab.

Gvalue gibt die Präzision an, mit der intern gearbeitet wird (bei keiner Angabe beträgt diese 400). Mittels der Mathematica Funktionen N und SetPrecision werden, falls nötig, die Eingaben der Dreieckskoordinaten auf die gewünschte Präzision gebracht.

Ein Test kann optional einen Hinweis darauf geben, wie vielen Stellen im Ergebnis wir vertrauen können. Dabei berechnen wir zusätzlich die Gewichte (OverTilde[w]) mit den Fehlerschranken (Subscript[ϵ, j]^2), j=1,2,3. Der Test terminiert, falls (Subscript[E, ϵ])=((||w-(OverTilde[w])(Subscript[||, 2]))/(||(OverTilde[w])(Subscript[||, 2]))) kleiner oder gleich einem in epstestvalue gewählten oder dem vorgebenen Wert (10^(-25)) ist. Für den Default-Wert (10^(-25)) liegen in der Datei epsMatrix annähernd geeignete Werte für (Subscript[ϵ, j]) vor. Wir erhalten diese Werte, indem wir die Quadraturen über die sphärischen Dreiecke mit den Ecken (Subscript[V, 1])=((Subscript[ρ, j]),0), (Subscript[V, 2])=((Subscript[ρ, j]),(Subscript[θ, k])), (Subscript[V, 3])=(0,0) berechnen, wobei (Subscript[ρ, j])=(10^(-j/8)), j=0,1,...,40, und (Subscript[θ, k])=(10^(-5))+(k/1000)π, k=0,1,...,49, (Subscript[θ, k])=(10^(-5))+((10(k-50)+50)/1000)π, k=50,51,...,140,  (Subscript[θ, k])=(10^(-5))+((k+9(k-141)+810)/1000)π, k=141,142,...,145, ist. Die Quadraturformeln sind dabei (näherungsweise) exakt für Spherical Harmonics  vom Höchstgrad n, 1≤ n ≤ 30.

Der Test bricht vorzeitig ab, falls ein oder mehrere Teildreiecke nicht berechnet werden (steuerbar über eps0value) oder der Algorithmus keine brauchbaren Werte für w=((Subscript[w, r,s])) liefert. Letzteres ist z.B. dann der Fall, wenn Gvalue zu gering ist, so dass zu viele Stellen ausgelöscht werden, oder die Werte  (Subscript[ϵ, 1]), (Subscript[ϵ, 2]), (Subscript[ϵ, 3 ]) zu groß sind. Der Quadratur können wir vertrauen, falls der Test einen Fehler (Subscript[E, ϵ])≠ -1 liefert, die Präzision p der Gewichte genügend groß und die Summe der Gewichte S=(∑_(k=1)^N)(Subscript[(w), k])≠0 ist (N sei die Anzahl der Gewichte und Knoten der Quadratur). Es kommt aber vor, dass der Gautschi-Algorithmus eine gewisse Gauss-Jacobi-Quadratur erfordert, deren Genauigkeitsgrad höher ist als die der vorberechneten Knoten im Verzeichnis Zeros bzw. Zeros_part. Es kann dann der Fall eintreten, dass bei einem höheren Polynomgrad, bei einem kleinen Radius und einem großen Winkel (Subscript[θ, k]) der Test zwar nicht abbricht, aber die berechneten Gewichte nicht mehr zuverlässig sind.  Die Absolutwerte der Gewichte werde überverhältnismäßig groß, die Genauigkeit der Summe S=(∑_(k=1)^N)(Subscript[(w), k]) nimmt ab.

Erste Schritte zur Berechnung einer Quadratur über einem sphärischen Dreieck:

Alle notwendigen Funktionen befinden sich im Package 

GaussianQuadratureSphericalTriangle.m, 

das an geeigneter Stelle der Mathematica Umgebung gespeichert sein sollte. Die Funktion GaussQuadratureST bestimmt näherungsweise die Gewichte w=((Subscript[w, r,s])) und Knoten x=((Subscript[χ, r]),(Subscript[ψ, s])), 1⩽|r|⩽n+2, 1⩽|s|⩽n+1, bei einem gegebenen sphärischen Dreieck Δ((Subscript[V, 1]),(Subscript[V, 2]),(Subscript[V, 3]))  und Polynomgrad n, so dass
 
    (∫_(Δ((Subscript[V, 1]), (Subscript[V, 2]), (Subscript[V, 3]))))(Subscript[Y, n])(x)  dξ(x)=(∑_(1⩽|r|⩽n+2))(∑_(1⩽|s|⩽n+1))(Subscript[w, r,s])(Subscript[Y, n])((Subscript[χ, r]),(Subscript[ψ, s])) für alle (Subscript[Y, n])∈(Subscript[ℍ, n]^2).

Mit der Funktion PlotTriangleNodes lassen sich das Dreieck und die Knoten plotten.
 
Zusätzlich benötigen wir ein Verzeichnis mit vorberechneten Nullstellen von orthogonalen Jacobi-Polynomen bzgl. der Gewichtsfunktionen (1-x)^r, r=-1/2 oder r=1/2 auf dem Intervall (-1,1). 
Angepasst an eine Arbeitspräzision von 400 und einen max. Polynomgrad n=30 für die zu bestimmende Quadratur liegen diese bis zum Polynomgrad Nmax=981 (r=-1/2) bzw. 
Nmaxtilde=648 (r=1/2) im Verzeichnis Zeros vorberechnet vor (Dateigröße etwa 153,82 MB); bis Polynomgrad 300 bzw. 100 im Verzeichnis Zeros_part (Dateigröße etwa 15,64 MB), dies reicht für die Berechnung der Beispiele in diesem Notebook aus. 
Im Notebook calculate_zeros.nb finden sich Funktionen, mit denen diese Nullstellen 
bestimmt werden können.

Als drittes (nach dem Package und einem Verzeichnis mit den Nullstellen) ist die Datei epsMatrix für unsere Beispielrechnungen notwendig. In ihr sind vorberechnete (Subscript[ϵ, j])-Werte (Fehlerschranken für das Abbruchkriterium im Gautschi-Algorithmus zur Berechnung von Dreiterm-Rekursions-Koeffizienten) gespeichert. Die vorberechneten Werte sind so gewählt, dass bei einem nicht zu hohen Präzisionsverlust die Quadraturgewichte auf mindestens 25 Nachkommastellen genau sind, bei einem Polynomgrad n ⩽ 30 als Maß der Exaktheit der Quadratur.

Wir nehmen an, dass sich das Package, das Verzeichnis Zeros_part und die Datei epsMatrix im aktuellen Mathematica Arbeitsverzeichnis oder einem Unterverzeichnis davon befinden.
1. Suchen des Pfades und Laden des Packages
2. Angabe des Pfades zu den vorberechneten Nullstellen.
3. Angabe der maximalen Nullstellen Nmax bzw. Nmaxtilde in Verzeichnis Zeros_part
4. Angabe des Pfades zu vorberechneten (Subscript[ϵ, j])-Werten in Datei epsMatrix

(*1.*)
myPathtoPackage=
FileNames["GaussianQuadratureSphericalTriangle.m",Directory[],Infinity][[1]];
pathtopackage=FileNameDrop[myPathtoPackage,-1];
Get["GaussianQuadratureSphericalTriangle`",Path->{pathtopackage}];

(*2.*)
pathzeros=FileNames["Zeros_part",Directory[],Infinity][[1]];

(*3.*)
Nmax=300;
Nmaxtilde=100;

(*4.*)
patheps=FileNames["epsMatrix",Directory[],Infinity][[1]];

Funktion GaussQuadratureST (Parameter, Optionen, Rückgabewerte) und 
Funktion PlotTriangleNodes (Plot des Dreiecks und der Knoten (optional))

GaussQuadratureST

Parameter: 
-Polynomgrad n
-die drei Ecken (Subscript[V, 1]),(Subscript[V, 2]),(Subscript[V, 3]) (gegen den Uhrzeigersinn angeordnet) des Dreiecks;
  in Form von Kugelkoordinaten (Subscript[V, j])=((Subscript[θ, j]),(Subscript[ρ, j])), j=1,2,3, oder
  in Form von kartesischen Koordinaten (Subscript[V, j])=((Subscript[v, j,1]),(Subscript[v, j,2]), (Subscript[v, j,3]))=(sin (Subscript[θ, j]) cos (Subscript[ρ, j]),sin (Subscript[θ, j]) sin (Subscript[ρ, j]),cos (Subscript[θ, j]))
 
 Optionen:
  Gvalue: gibt die Arbeitspräzision an, d.h. intern wird mit dieser Präzision gerechnet.
                Default-Wert ist 400 (bei Angabe eines höheren Wertes können unter                      
                Umständen Probleme auftauchen)
 eps0value: Es sei γ der Winkel  gegenüber der Basis in einem gleichschenkligen Dreieck. 
                    Ist |π-γ| ≤ 10*eps0value, dann wird für dieses Dreieck die Quadratur nicht berechnet.
                    Default-Wert ist (10^(-16)).
funct: Name des Algorithmus zur Berechnung von Dreiterm-Rekursions-Koeffizienten. Es steht ein
           diskretisierter Stieltjes-Algorithmus oder ein Lanczos-Algorithmus zur Wahl. 
           Default-Wert ist "lanczos" (alternativ kann "stieltjes" gewählt werden).
 epslistvalue: ist eine Liste mit den Einträgen (Subscript[ϵ, 1]), (Subscript[ϵ, 2]), (Subscript[ϵ, 3]). Wird ein sphärisches Dreieck in drei
                        gleichschenklige Dreiecke (Subscript[Δ, j]),j=1,2,3, aufgeteilt (z.B. weil es selber nicht 
                        gleichschenklig ist), dann wird für jedes dieser Teildreiecke die Quadratur bestimmt. 
                        Mit dem Gautschi-Algorithmus werden dabei Approximationen von Dreiterm-
                        Rekursions-Koeffizienten berechnet. Für das jeweilige  Abbruchkriterium benötigen
                        wir eine Fehlerschrank (Subscript[ϵ, j]) (erfolgt keine Aufteilung in Teildreiecke, dann ist nur (Subscript[ϵ, 1]) 
                        relevant).
                        Default-Liste ist {}; es wird aus vorberechneten Werten in der Datei epsMatrix ein
                        annähernd geeigneter Wert genommen (besonders geeignet bei epstestvalue-
                        Default-Wert (10^(-25)) und 1 ≤ n ≤ 30).
 epstestboolean: gibt an, ob ein Test durchgeführt werden soll. 
                              Default-Wert ist False (epstestboolean->True: Test wird durchgeführt).
 epstestvalue:   gibt die Fehlerschranke im Test an. Falls der Fehler (Subscript[E, ϵ])=((||w-(OverTilde[w])(Subscript[||, 2]))/(||(OverTilde[w])(Subscript[||, 2])))  sie
                            unterschreitet, wird der Testlauf abgebrochen. Die Einträge in (OverTilde[w]) werden
                            mit den Fehlerschranken (Subscript[ϵ, j]^2), j=1,2,3, im Gautschi-Algorithmus berechnet.
                            Default-Wert ist (10^(-25)).
 partboolean:   gibt an, ob ein gleichschenkliges Dreieck (wie andere Dreiecke auch) in 
                          wiederum  gleichschenklige Dreiecke aufgeteilt wird. 
                          Default-Wert ist False (epstestboolean->True: Dreieck wird aufgeteilt).                         
epsPivalue: Es sei ρ der Radius des Dreiecksumkreises. Ist |π/2-ρ| ≤ epsPivalue, dann tritt 
                     der Spezialfall ρ=π/2 ein.  
                     Default-Wert: (10^(-16))                       
    
Rückgabewerte:
-Liste der Näherungen der Gewichte w=((Subscript[w, k])), k=1,...,4(Subscript[n, 4])(Subscript[n, 5]), (Subscript[n, 4])=n+1, (Subscript[n, 5])=n+2
-Liste der Näherungen der Knoten x=( (Subscript[x, k]))=(((Subscript[ψ, k]),(Subscript[χ, k]))), k=1,...,4(Subscript[n, 4])(Subscript[n, 5]), (Subscript[n, 4])=n+1, (Subscript[n, 5])=n+2, mit den
  äquatorialen Winkeln (Subscript[ψ, k]), 0⩽(Subscript[ψ, k])<π  und den polaren Winkeln (Subscript[χ, k]), 0⩽(Subscript[χ, k])⩽π
-Zeit t der Berechnung (ohne Test) 
-Fehlerschranken (Subscript[ϵ, 1]), (Subscript[ϵ, 2]), (Subscript[ϵ, 3]) im Abbruchkriterium des Gautschi-Algorithmus, mit denen die
 Fehler der Quadraturformeln über drei Teildreiecke (Subscript[Δ, j]) beeinflusst werden kann, in das ein bel.
  Dreieck gegebenfalls aufgeteilt wird (wird es nicht aufgeteilt, dann ist (Subscript[ϵ, 1])=(Subscript[ϵ, 2])=(Subscript[ϵ, 3]))
-Liste (Subscript[N, L])mit der Anzahl der Gewichte und Knoten der Gauss-Jacobi-Quadraturen zur Approximation
  von Dreiterm-Rekursions-Koeffizienten im Gautschi-Algorithmus.
-Werte (Subscript[n, 4]) und (Subscript[n, 5]): (Subscript[n, 5])=n+2, (Subscript[n, 4])=n+1; Anzahl der Gewichte und Knoten der Quadratur N=4(Subscript[n, 4])(Subscript[n, 5])
-Flächeninhalt F  des Dreiecks 
-Fehler (Subscript[E, ϵ])=((||w-(OverTilde[w])(Subscript[||, 2]))/(||(OverTilde[w])(Subscript[||, 2]))) berechnet bei Durchlauf des Tests; Einträge in (OverTilde[w]) werden mit Fehlerschranken
   (Subscript[ϵ, j]^2), j=1,2,3, im Gautschi-Algorithmus berechnet. Wird kein Test durchgeführt oder bricht Test
  vorzeitig ab, dann ist   (Subscript[E, ϵ])=-1.
-Anzahl Z der Testdurchläufe 


      

PlotTriangleNodes:
  
Parameter:
 die drei Ecken (Subscript[V, 1]), (Subscript[V, 2]), (Subscript[V, 3]) (gegen den Uhrzeigersinn angeordnet) des Dreiecks; 
 in Form von Kugelkoordinaten (Subscript[V, j])=((Subscript[θ, j]),(Subscript[ρ, j])), j=1,2,3, oder
 in Form von kartesischen Koordinaten (Subscript[V, j])=((Subscript[v, j,1]),(Subscript[v, j,2]),(Subscript[v, j,3]))=(sin (Subscript[θ, j]) cos (Subscript[ρ, j]),sin (Subscript[θ, j]) sin (Subscript[ρ, j]),cos (Subscript[θ, j]))

Optionen:
 Nodes: gibt eine Liste mit Punkten in polaren Koordinaten an, die geplottet werden sollen.
              Default-Matrix ist {{}}.
 labelboolean: gibt an, ob die Ecken, der Nordpol und der Mittelpunkt des Umkreises mit   
                         Bezeichnungen versehen werden und der Umkreis geplottet wird. 
                         Default-Wert ist True  
                         (labelboolean->False: keine Bezeichnungen, kein Plot des Umkreises).

Returns : eine dreidimensionale Graphik der Sphäre mit einem sphärischen Dreieck und
                geplotteten Punkten (falls gewünscht)

1. Beispiel:
Dreieck mit Plot (Ecken (Subscript[V, 1]),(Subscript[V, 2]),(Subscript[V, 3]), Nordpol (Subscript[n, 0]), Umkreis des Dreiecks mit Mittelpunkt m)

n=4;
(Subscript[θ, 1])=Pi/3;
(Subscript[ρ, 1])=Pi/10+Pi;
(Subscript[θ, 2])=3*Pi/4;
(Subscript[ρ, 2])=Pi/7+Pi;
(Subscript[θ, 3])=Pi/4;
(Subscript[ρ, 3])=Pi/5+Pi;


(Subscript[t, 2])=AbsoluteTiming[{w,x,t,(Subscript[ϵ, 1]),(Subscript[ϵ, 2]),(Subscript[ϵ, 3]),(Subscript[N, L]),(Subscript[n, 5]),(Subscript[n, 4]),F,(Subscript[E, ϵ]),Z}=GaussQuadratureST[n,{(Subscript[θ, 1]),(Subscript[ρ, 1])},{(Subscript[θ, 2]),(Subscript[ρ, 2])},{(Subscript[θ, 3]),(Subscript[ρ, 3])}];];

Print["Gewichte:"];
N[w]

PlotTriangleNodes[{(Subscript[θ, 1]),(Subscript[ρ, 1])},{(Subscript[θ, 2]),(Subscript[ρ, 2])},{(Subscript[θ, 3]),(Subscript[ρ, 3])},Nodes->x]

Zusätzliche Informationen zum ersten Beispiel:

N=Anzahl der Gewichte und Knoten: N=4(Subscript[n, 4])(Subscript[n, 5])
F=Flächeninhalt des Dreiecks
S=Summe der Gewichte: S=(∑_(k=1)^N)(Subscript[w, k])
Fehler (Subscript[E, F=])((|S-F|)/F)
p=Präzision der Gewichte w
t=Zeit für die Berechnung der Quadraturformel ohne Test und Vorberechnungen in sec
(Subscript[t, 2])=Gesamtberechnungszeit für die Berechnung der Quadraturformel mit Vorberechnungen (und 
     evt. mit Test) in sec
(Subscript[ϵ, j])=Werte, die die Genauigkeit der Quadraturformel der Teildreiecke beeinflussen (j=1,2,3)
Z=Anzahl der Testdurchläufe
Fehler (Subscript[E, ϵ])=((||w-(OverTilde[w])(Subscript[||, 2]))/(||(OverTilde[w])(Subscript[||, 2]))) im Test. Die Einträge in (OverTilde[w]) werden mit Fehlerschranken (Subscript[ϵ, j]^2), j=1,2,3,  im 
                              Abbruchkriterium des Gautschi-Algorithmus berechnet. Wird kein Test 
                              durchgeführt oder bricht der Test  vorzeitig ab, dann ist (Subscript[E, ϵ])=-1.
(Subscript[N, L,max=])= Maximum der Anzahl der Gewichte und Knoten der Gauss-Jacobi-Quadraturen zur        
               Approximation von Dreiterm-Rekursions-Koeffizienten im Gautschi-Algorithmus

S=Plus@@w;
(Subscript[E, F])=Abs[S-F]/F;
TableForm[{4*(Subscript[n, 4])*(Subscript[n, 5]),N[F,25],N[S,25],N[(Subscript[E, F]),25],Precision[w],t},TableHeadings->{{"N","F","S","(Subscript[E, F])","p","t/sec"}}]

TableForm[{(Subscript[t, 2])[[1]],-Log[10,(Subscript[ϵ, 1])],-Log[10,(Subscript[ϵ, 2])],-Log[10,(Subscript[ϵ, 3])],N[(Subscript[E, ϵ]),25],Z,Max[(Subscript[N, L])]},TableHeadings->
{{"(Subscript[t, 2])/sec","-(Subscript[log, 10])(Subscript[ϵ, 1])","-(Subscript[log, 10])(Subscript[ϵ, 2])","-(Subscript[log, 10])(Subscript[ϵ, 3])","(Subscript[E, ϵ])","Z","(Subscript[N, L, max])"}}]

Hier wird die Berechnung wiederholt, aber wir wählen die Arbeitspräzision Gvalue=50 und die Fehlerschranke $MachineEpsilon in den Optionen epslistvalue und epstestvalue. Option epstestboolean wird auf True gesetzt, so dass der Test ausgeführt wird. Die Berechnung liefert:

N=Anzahl der Gewichte und Knoten: N=4(Subscript[n, 4])(Subscript[n, 5])
F=Flächeninhalt des Dreiecks
S=Summe der Gewichte: S=(∑_(k=1)^N)(Subscript[w, k])
Fehler (Subscript[E, F=])((|S-F|)/F)
p=Präzision der Gewichte w
t=Zeit für die Berechnung der Quadraturformel ohne Test und Vorberechnungen in sec
(Subscript[t, 2])=Gesamtberechnungszeit für die Berechnung der Quadraturformel mit Vorberechnungen (und 
     evt. mit Test) in sec
(Subscript[ϵ, j])=Werte, die die Genauigkeit der Quadraturformel der Teildreiecke beeinflussen (j=1,2,3)
Z=Anzahl der Testdurchläufe
Fehler (Subscript[E, ϵ])=((||w-(OverTilde[w])(Subscript[||, 2]))/(||(OverTilde[w])(Subscript[||, 2]))) im Test. Die Einträge in (OverTilde[w]) werden mit Fehlerschranken (Subscript[ϵ, j]^2), j=1,2,3,  im 
                              Abbruchkriterium des Gautschi-Algorithmus berechnet. Wird kein Test 
                              durchgeführt oder bricht der Test  vorzeitig ab, dann ist (Subscript[E, ϵ])=-1.
(Subscript[N, L,max=])= Maximum der Anzahl der Gewichte und Knoten der Gauss-Jacobi-Quadraturen zur        
               Approximation von Dreiterm-Rekursions-Koeffizienten im Gautschi-Algorithmus

(Subscript[t, 2])=AbsoluteTiming[{w,x,t,(Subscript[ϵ, 1]),(Subscript[ϵ, 2]),(Subscript[ϵ, 3]),(Subscript[N, L]),(Subscript[n, 5]),(Subscript[n, 4]),F,(Subscript[E, ϵ]),Z}=GaussQuadratureST[n,{(Subscript[θ, 1]),(Subscript[ρ, 1])},{(Subscript[θ, 2]),(Subscript[ρ, 2])},{(Subscript[θ, 3]),(Subscript[ρ, 3])},Gvalue->50,epslistvalue->{$MachineEpsilon,$MachineEpsilon,$MachineEpsilon},epstestvalue->$MachineEpsilon,epstestboolean->True];];



S=Plus@@w;

(Subscript[E, F])=Abs[S-F]/F;
TableForm[{4*(Subscript[n, 4])*(Subscript[n, 5]),N[F],N[S],N[(Subscript[E, F])],Precision[w],t},TableHeadings->{{"N","F","S","(Subscript[E, F])","p","t/sec"}}]


TableForm[{(Subscript[t, 2])[[1]],-Log[10,(Subscript[ϵ, 1])],-Log[10,(Subscript[ϵ, 2])],-Log[10,(Subscript[ϵ, 3])],(Subscript[E, ϵ]),Z,Max[(Subscript[N, L])]},TableHeadings->{{"(Subscript[t, 2])/sec","-(Subscript[log, 10])(Subscript[ϵ, 1])","-(Subscript[log, 10])(Subscript[ϵ, 2])","-(Subscript[log, 10])(Subscript[ϵ, 3])","(Subscript[E, ϵ])","Z","(Subscript[N, L, max])"}}]

2. Beispiel:
Berechnung der Quadratur und Plot des Dreiecks und der Knoten
(Ecken (Subscript[V, 1]),(Subscript[V, 2]),(Subscript[V, 3]), Nordpol (Subscript[n, 0]), Dreicksumkreis mit Zentrum m)

n=8;
(Subscript[θ, 1])=Pi/3;
(Subscript[ρ, 1])=Pi/4;
(Subscript[θ, 2])=3*Pi/4;
(Subscript[ρ, 2])=3*Pi/4;
(Subscript[θ, 3])=Pi/4;
(Subscript[ρ, 3])=Pi/5+2*Pi/3;


(Subscript[t, 2])=AbsoluteTiming[{w,x,t,(Subscript[ϵ, 1]),(Subscript[ϵ, 2]),(Subscript[ϵ, 3]),(Subscript[N, L]),(Subscript[n, 5]),(Subscript[n, 4]),F,(Subscript[E, ϵ]),Z}=GaussQuadratureST[n,{(Subscript[θ, 1]),(Subscript[ρ, 1])},{(Subscript[θ, 2]),(Subscript[ρ, 2])},{(Subscript[θ, 3]),(Subscript[ρ, 3])}];];
Print["Gewichte:"];
N[w]

PlotTriangleNodes[{(Subscript[θ, 1]),(Subscript[ρ, 1])},{(Subscript[θ, 2]),(Subscript[ρ, 2])},{(Subscript[θ, 3]),(Subscript[ρ, 3])},Nodes->x]

S=Plus@@w;
(Subscript[E, F])=Abs[S-F]/F;
TableForm[{4*(Subscript[n, 4])*(Subscript[n, 5]),N[F,25],N[S,25],N[(Subscript[E, F]),25],Precision[w],t},TableHeadings->{{"N","F","S","(Subscript[E, F])","p","t/sec"}}]



TableForm[{(Subscript[t, 2])[[1]],-Log[10,(Subscript[ϵ, 1])],-Log[10,(Subscript[ϵ, 2])],-Log[10,(Subscript[ϵ, 3])],N[(Subscript[E, ϵ]),20],Z,Max[(Subscript[N, L])]},TableHeadings->{{"(Subscript[t, 2])/sec","-(Subscript[log, 10])(Subscript[ϵ, 1])","-(Subscript[log, 10])(Subscript[ϵ, 2])","-(Subscript[log, 10])(Subscript[ϵ, 3])","(Subscript[E, ϵ])","Z","(Subscript[N, L, max])"}}]

3. Beispiel
Oktant der Sphäre und Plot der Oktanten mit den Quadraturknoten

n=4;
X0={0,0,1};
alpha=0;
beta=Pi/2;
P1={Cos[alpha],Sin[alpha],0};
P2={Cos[beta],Sin[beta],0};

(Subscript[t, 2])=AbsoluteTiming[{w,x,t,(Subscript[ϵ, 1]),(Subscript[ϵ, 2]),(Subscript[ϵ, 3]),(Subscript[N, L]),(Subscript[n, 5]),(Subscript[n, 4]),F,(Subscript[E, ϵ]),Z}=GaussQuadratureST[n,X0,P1,P2,partboolean->True];];


S=Plus@@w;
(Subscript[E, F])=Abs[S-F]/F;
TableForm[{4*(Subscript[n, 4])*(Subscript[n, 5]),N[F,25],N[S,25],N[(Subscript[E, F]),25],Precision[w],t},TableHeadings->{{"N","F","S","(Subscript[E, F])","p","t/sec"}}]



TableForm[{(Subscript[t, 2])[[1]],-Log[10,(Subscript[ϵ, 1])],-Log[10,(Subscript[ϵ, 2])],-Log[10,(Subscript[ϵ, 3])],N[(Subscript[E, ϵ]),20],Z,Max[(Subscript[N, L])]},TableHeadings->{{"(Subscript[t, 2])/sec","-(Subscript[log, 10])(Subscript[ϵ, 1])","-(Subscript[log, 10])(Subscript[ϵ, 2])","-(Subscript[log, 10])(Subscript[ϵ, 3])","(Subscript[E, ϵ])","Z","(Subscript[N, L, max])"}}]

Show[(xtemp=N[x,25];
xtemp[[All,1]]=N[x[[All,1]],25]+#*Pi/2;
alpha=#*Pi/2;
beta=Pi/2+#*Pi/2;
P1={Cos[alpha],Sin[alpha],0};
P2={Cos[beta],Sin[beta],0};
plot1=PlotTriangleNodes[X0,P1,P2,Nodes->xtemp,labelboolean->False];
xtemp[[All,2]]=Pi-N[x[[All,2]],25];
plot2=PlotTriangleNodes[-X0,P2,P1,Nodes->xtemp,labelboolean->False];
Show[plot1,plot2])&/@{0,1,2,3}]

Testen,ob die Quadratur exakt ist für Spherical Harmonics vom Höchstgrad n,d.h. ob
(∫_(Δ((Subscript[V, 1]), (Subscript[V, 2]), (Subscript[V, 3]))))(Subscript[Y, n])(x)  dξ(x)=(∑_(1⩽|r|⩽n+2)) (∑_(1⩽|s|⩽n+1))(Subscript[w, r,s])(Subscript[Y, n])((Subscript[χ, r]),(Subscript[ψ, s]))=(Sqrt[4π])(Subscript[δ, n,0])füralle(Subscript[Y, n])∈(Subscript[ℍ, n]^2):

TestQuadratureOctant[n_,G_,WeigthsList_,NodesList_]:=
Module[{GMatrix,Ytemp,xcoord,ycoord,ktemp,bvector,expvektor,summe,m,l},
GMatrix=ConstantArray[0,(n+1)^2];GMatrix[[1]]=Table[N[1/(2*Sqrt[Pi]),G],{j,1,Length[NodesList[[All,1]]]}];

(Ytemp=(-1)^(#[[2]])*SphericalHarmonicY[#[[1]],#[[2]],xcoord,ycoord]/.{xcoord->N[NodesList[[All,2]],G],ycoord->N[NodesList[[All,1]],G]};
GMatrix[[(#[[1]]+1)^2-(#[[1]]-#[[2]])]]=Ytemp;GMatrix[[(#[[1]]+1)^2-(#[[1]]+#[[2]])]]=Conjugate[Ytemp])&/@Flatten[Table[{l,m},{l,1,n},{m,0,l}],1];


bvector=GMatrix.WeigthsList;
expvektor=Exp[I*#*ktemp*Pi/2]& /@Flatten[Table[m,{l,0,n},{m,-l,l}]];


summe=Plus@@(
(IdentityMatrix[(n+1)^2]+(DiagonalMatrix[If[Mod[#[[1]]-Abs[#[[2]]],2]==0,1,-1]&/@Flatten[Table[{l,m},{l,0,n},{m,-l,l}],1]])).((expvektor/.{ktemp->#})*bvector)&/@{0,1,2,3});


Return[summe];

];

MatrixForm[TestQuadratureOctant[n,25,w,x]]

Literaturverzeichnis

[1]  J. Beckmann, H. N. Mhaskar, J. Prestin, Quadrature formulas for integration of multivariate   
       trigonometric polynomials on spherical triangles, GEM Int. J. Geomath. 3(1) (2012), 119-138
 [2]  J. Beckmann, H. N. Mhaskar, J. Prestin, Local numerical integration on the sphere, GEM 
       Int. J. Geomath. 5(2) (2014), 143-162
 [3]  W. Gautschi, On Generating Orthogonal Polynomials, SIAM J. Sci. Stat. Comput. 3(3) (1982), 
       289-317
 [4]  W. Gautschi, Orthogonal Polynomials: Computation and Approximation, Oxford University 
       Press, Oxford, 2004
 [5]  W. Gautschi, Orthogonal polynomials (in Matlab), J. Comput. Appl. Math. 178(1-2) (2005), 
       215-234  
 [6]  W.  Gautschi, Package OPQ: a Matlab Suite of Programs for Generating Orthogonal 
       Polynomials and Related Quadrature Rules,
       http:// www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html
 [7]  Q. T. Le Gia, H. N. Mhaskar, Localized linear polynomial oprerators and quadrature formulas
       on the sphere, SIAM J. Numer. Anal. 47(1) (2008), 440-466
 [8]  H. N. Mhaskar, F. J. Narcowich, J. D. Ward, Spherical Marcinkiewcz-Zygmund inequalities and  
       positive quadrature, Math. Comp. 70(235) (2001), 1113-1130     
 [9]  H. N. Mhaskar, Local quadrature formulas on the sphere, J. Complexity 20(5) (2004), 753-772  
[10]  H. N. Mhaskar, Local quadrature formulas on the sphere, II, in “Advances in Constructive 
       Approximation” (M. Neamtu and E. B. Saff eds.), Nashboro Press, Nashville, 2004, pp.333-344
[11] R. J. Renka, Algortihm 772:STRIPACK:Delaunay triangulation and Voronoi diagram on the
       surface of a sphere, ACM Trans. Math. Softw. 23(3) (1997), 416-434