ガウス積分について

質問

現在数値計算を行っているのですが、積分点数をかなり増やす必要がでてきました。
最低10個程度は必要なのでn個の積分点の位置と重みを計算するCプログラムを探しているのですが、見つかりません。
検索して出てくるのは、大抵5個までで、プログラム内に直接数値を打ち込んでいるものでした。
積分点の個数を与えると、各位置と重みを出力するプログラムを教えてください。

前回の質問はどうやらカテゴリーが不適だったようなので、再度質問させて頂きました。

回答No.2

f272 · ベストアンサー

http://www7.ocn.ne.jp/~kawa1/numeric.pdf
の204ページの公式(5.112)をエクセルVBAで書いた例

Sub gg()
Dim XW() As Double
Dim i As Long, j As Long, m As Long
Dim p1 As Double, p2 As Double, p3 As Double, pp As Double
Dim z As Double, z1 As Double

n = Range("a1") '積分点の個数
ReDim XW(1 To n, 1 To 2) As Double

m = (n + 1) / 2
For i = 1 To m
z = Cos(3.141592654 * (i - 0.25) / (n + 0.5))
Do
p1 = 1
p2 = 0
For j = 1 To n
p3 = p2
p2 = p1
p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j
Next j
pp = n * (z * p1 - p2) / (z * z - 1)
z1 = z
z = z1 - p1 / pp
Loop While (Abs(z - z1) > 0.00000000001)
XW(i, 1) = z
XW(n + 1 - i, 1) = z
XW(i, 2) = 2 / ((1 - z * z) * pp * pp)
XW(n + 1 - i, 2) = XW(i, 2)
Next i

Range("a2").Resize(n, 2) = XW '各位置(A列)と重み(B列)

End Sub

これをCにするくらいはどうってことないよね。

お礼

do～loop whileで苦労しましたが、なんとかCに変換でき、
積分点・重みを求めることができました。
回答ありがとうございました。

回答No.1

FT56F001

プログラムが欲しいという回答になっていませんが，数学板で拝見したので回答します。
ガウスの分点法で数値計算したいというご要望ですね。

有限区間を数値積分する，ガウスの分点法の分点と重みは，
ルジャンドル多項式のゼロ点とその微分で計算できます。
具体的には，計算する積分区間を[-1,1]に変換した後，
n次のルジャンドル多項式P(x)のゼロ点をx_i，ルジャンドル多項式の微分をP'(x)として，
A_i=(1/P'(x))∫[x=-1から1]P(x)/(x-x_i)dx
で重みA_iを計算できます。これで求まるn個の分点と重みで，
2n-1次以下の多項式を誤差なく積分できます。

ただし，どんな積分に利用されるのかわかりませんが，関数値を1点計算するのに，
とても手間がかかるケースなのでしょうか?
一般論として言うと，高次のガウスの方法のような凝った積分公式よりも，
単純な中点公式か台形公式，せいぜいシンプソンの公式で，
単純に分点を増やして細かく計算する方が，
プログラムが簡単で，たくましい計算法になります。
(たくましいとは，多項式のような素直な関数だけでなく，
　区分連続など条件が悪い関数でも，それなりに積分できるという意味です)

昨今，計算機がパワフルなので，人が凝ったプログラムを作る時間より，
単純な積分法で細かく分割して計算機をぶん回す方が，結局は早いのではないかと思っています。

お礼

おいおいは変更することもあるかもしれませんが、
現在はガウス積分でいくようです。
あまりにも多くの積分点が必要になるようなら変更するかもしれません。
回答ありがとうございました。