Automatische Erzeugung von Integrationspunkten und Gewichten für Dreiecke und Tetraeder

12

Normalerweise konsultiert man ein Papier oder ein Buch, um Integrationspunkte und -gewichte für Einheitsdreieck und Tetraeder zu finden. Ich suche nach einer Methode, um solche Punkte und Gewichte automatisch zu berechnen. Das folgende Mathematica- Codebeispiel berechnet Integrationsgewichte und -punkte für Einheitenlinienelemente (Quad / Hexaeder):

unitGaussianQuadraturePoints[points_] := 
  Sort[x /. 
    Solve[Evaluate[LegendreP[points, x] == 0], {x}], ! 
     OrderedQ[N[{#1, #2}]] &];

unitGaussianQuadratureWeights[points_] := 
  Module[{gps, f, int, integr, vars, eqns}, 
   gps = unitGaussianQuadraturePoints[points];
   f[0, 0] := 1;
   f[0., 0] := 1.;
   f[x_, n_] := x^n;
   int = Integrate[f[x, #], x] & /@ Range[0, points - 1];
   integr = Subtract @@@ (int /. x :> {1, -1});
   vars = Table[Unique[c], {Length[gps]}];
   eqns = 
    Table[Plus @@ Thread[Times[vars, f[#, i - 1] & /@ gps]] == 
      integr[[i]], {i, points}];
   Return[(vars /. Solve[eqns, vars])];];


unitGaussianQuadratureWeights[2]

{{1, 1}}

unitGaussianQuadraturePoints[2]

{1/Sqrt[3], -(1/Sqrt[3])}

Ich suche eine Arbeit / ein Buch, das algorithmisch beschreibt, wie dies für Dreiecke und / oder Tetraeder gemacht wird. Kann mich jemand auf einige Informationen dazu hinweisen. Vielen Dank.

David Ketcheson
quelle
1
Es gibt einen einfacheren Weg , um Ihre Gauß-Legendre Quadratur Regeln zu tun Mathematica : {points, weights} = MapThread[Map, {{2 # - 1 &, 2 # &}, Most[NIntegrate`GaussRuleData[n, prec]]}].
JM
Auf jeden Fall: Sie gesehen haben diese ?
JM
@JM, Ihre oben vorgeschlagene Methode funktioniert leider nicht für prec = Infinity; aber danke auch dafür.
2
In diesem Fall funktioniert die folgende Methode dank Golub und Welsch: Transpose[MapAt[2(First /@ #)^2 &, Eigensystem[SparseArray[{Band[{2, 1}] -> #, Band[{1, 2}] -> #}, {n, n}]], {2}]] &[Table[k/Sqrt[(2 k - 1)(2 k + 1)], {k, n - 1}]] .
JM
1
Hier ist das Papier von Golub und Welsch. Ich werde meine Papiere durchforsten und sehen, ob es etwas für Vereinfachungen gibt ...
JM

Antworten:

3

Hier ist ein Papier Artikel http://journal.library.iisc.ernet.in/vol200405/paper6/rathod.pdf , in dem beschrieben wird, wie das Einheitendreieck auf das Standard-2-Quadrat abgebildet wird, um die Gewichte und Stichprobenpunkte für die Dreieck in Form von Gauß-Legendre-Punkten für das Standard-2-Quadrat.

James Custer
quelle
Das ist eine interessante Idee, es sieht so aus, als ob für n = 2 4 Punkte benötigt werden, für die typische Literaturreferenz für Dreiecke für n = 2 werden 3 Punkte angegeben. Wissen Sie etwas darüber?
Dies ist darauf zurückzuführen, dass sie eine Zuordnung vom Dreieck zum Quadrat verwenden. Darüber hinaus kann ich nichts sagen, da ich nicht mit Dreiecken arbeite (ich benutze Vierecke), also weiß ich nicht, was in der Praxis normalerweise gemacht wird. Ich habe gerade die Zeitung gefunden und dachte, es wäre ziemlich einfach.
James Custer
Es ist in der Tat ganz einfach und ich werde sehen, dass die anderen Artikel dies vorschlagen, aber die Einfachheit dieses Artikels und die Eleganz, etwas zu verwenden, das ich bereits habe, sind ein Plus für diesen Artikel. Der Nachteil ist dann die zusätzliche Funktionsbewertung. Danke auf jeden Fall.
Ein weiterer Nachteil ist, dass die Punkte nicht symmetrisch sind.