Automatyczne generowanie punktów integracji i ciężarów dla trójkątów i czworościanów

12

Zwykle należy skonsultować się z gazetą lub książką, aby znaleźć punkty integracji i ciężary dla trójkąta jednostkowego i czworościanów. Szukam metody automatycznego obliczania takich punktów i wag. Poniższy przykład kodu Mathematica oblicza wagi integracji i punkty dla elementów linii jednostkowej (quad / hexahedron):

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])}

Szukam papieru / książki, która opisuje algorytmicznie, jak to się robi w przypadku trójkątów i / lub czworościanów. Czy ktoś może wskazać mi jakieś informacje na ten temat. Dzięki.

David Ketcheson
źródło
1
Jest łatwiejszy sposób, aby robić swoje reguły kwadratury Gaussa-Legendre'a w Mathematica : {points, weights} = MapThread[Map, {{2 # - 1 &, 2 # &}, Most[NIntegrate`GaussRuleData[n, prec]]}].
JM
W każdym razie: widziałeś to ?
JM
@JM, twoja wyżej zaproponowana metoda niestety nie działa dla prec = Infinity; ale dzięki za to.
2
W takim przypadku, tutaj jest to metoda, która, ze względu na prace Golub i 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
Oto artykuł Goluba i Welscha. Przekopię się przez papiery i sprawdzę, czy jest coś prostego ...
JM

Odpowiedzi:

3

Oto papier http://journal.library.iisc.ernet.in/vol200405/paper6/rathod.pdf, który opisuje, w jaki sposób odwzorować trójkąt jednostkowy na standardowy 2-kwadratowy w celu obliczenia ciężarów i punktów próbkowania dla trójkąt pod względem punktów Gaussa-Legendre'a dla standardowego kwadratu 2.

James Custer
źródło
To ciekawy pomysł, wygląda na to, że dla n = 2 potrzeba 4 punktów, dla typowego odniesienia literaturowego dla trójkątów dla n = 2 podano 3 punkty. Wiesz cokolwiek o tym?
Wynika to z faktu, że używają odwzorowania z trójkąta na kwadrat. Nie mogę nic więcej powiedzieć, ponieważ nie pracuję z trójkątami (używam czworoboków), więc nie wiem, co zwykle robi się w praktyce. Właśnie znalazłem gazetę i pomyślałem, że to całkiem prosta rzecz do zrobienia.
James Custer
w rzeczy samej jest to dość proste i zobaczę, że sugerują to inne artykuły, ale prostota tego i elegancja korzystania z czegoś, co już mam, są zaletą tego. Minusem jest zatem ocena funkcji dodatkowej. Dzięki w każdym razie.
kolejną wadą jest to, że punkty nie są symetryczne.