プログラミング演習（Pint）数値積分

１次元の関数 f(x) の定積分$\int_a^b f(x) dx$ の数値積分のうち、台形公式とガウスルジャンドル積分について学ぶ。

(理論)台形公式

区間 [a,b] を N 当分し、それぞれの小区間の面積を台形で近似する。 h=(b-a)/N として、定積分は

$$\int_a^b f(x) dx \sim \frac{h}{2}[f(a)+f(a+h)]+\frac{h}{2}[f(a+h)+f(a+2h)]+ ... *\frac{h}{2}[f(a+(N-1)h)+f(b)]$$$$=\left[\frac{1}{2}\{f(a)+f(b)\}+f(a+h)+...f(a+(N-1)h)\righ$$

と近似される。





（例題Pintー１）下のプログラムは、台形公式により １０当分して x⁴　を 1 から 2まで積分するプログラムです。作成して実行しましょう。



gcc pexint-1.c -o pexint-1 -lm
./pexint-1
       6.2002333330000
analytical value 6.2000000000000

（課題Pint-１）
例題を参考にして、台形公式により f(x)=4/(1+x²) を a= 0 から b=1 まで積分するプログラムを作成し、実行しなさい。分割数を4, 8, 16, 32, 64, 128, 256, 512, ... としたとき、積分の結果の収束性を調べなさい。この積分は解析的に解ける。台形公式における結果と解析解とを比較せよ。

（課題Pint-２)
台形公式により f(x)=exp(-|x|) について a から b まで積分するプログラムを作成しなさい。a, b, と分割数 n については入力するようにして、また、f(x) について関数を用いること。この積分は解析的に解ける。台形公式における結果と解析解とを比較せよ。 a<0 かつ b>0 について試してみること。たとえば a=-1.111111111, b=1.0 など。

(理論)ガウス・ルジャンドル積分
関数f(x)をルジャンドル多項式展開で近似し、そのn次までの積分値で積分を近似する方法をガウス・ルジャンドル積分という。ルジャンドル多項式展開した式の積分を評価する数学的導出は省略するが、最終的に以下のような近似公式を得る：
$${\int }_a^{b}f(x)dx\simeq \sum _{i=1}^N{w}_{i}f(x_i)$$ ここで、w_iは重みと呼ばれる量であり、x_iはガウス点呼ばれる特定のxの値である。台形公式などは、重みw_iが均等で、x_iも等間隔に取っている積分公式と見なせる。ガウス・ルジャンドル積分は、重みw_iとガウス点x_iを非対等・非等間隔に上手くとって、台形公式よりも効率的に積分を評価する方法だと言える。

区間[a,b]とガウス点の数Nが与えられれば、(w_i,x_i)の値はfによらずに決まる。そのため、[a,b]とガウス点の数Nから(x_i,w_i)を計算するパッケージは、様々な数値計算ライブラリに実装されており、自分で(x_i,w_i)を計算するコードを書く必要は無い。例えば[0,1]区間の積分N個のガウス点で評価する際の(x_i,w_i)は以下のファイルのように与えられる(１列目がx_i, ２列目がw_i)：
GLwx_N1.txt (N=1の時のx_iとw_i)
GLwx_N2.txt (N=2の時のx_iとw_i)
GLwx_N4.txt (N=4の時のx_iとw_i)
GLwx_N8.txt (N=8の時のx_iとw_i)
GLwx_N16.txt (N=16の時のx_iとw_i)
GLwx_N32.txt (N=32の時のx_iとw_i)
GLwx_N64.txt (N=64の時のx_iとw_i)
GLwx_N128.txt (N=128の時のx_iとw_i)
GLwx_N256.txt (N=256の時のx_iとw_i)
GLwx_N512.txt (N=512の時のx_iとw_i)
GLwx_N1024.txt (N=1024の時のx_iとw_i)
ちなみに上記の和を計算する際に通常通り和をとっても良いし、 $$\overrightarrow{w}=[w_1,w_2,...,w_N]\overrightarrow{f}=[f(x_1),f(x_2),...,f(x_N)]$$ というベクトルを準備しておいて、 $${\displaystyle {\int }_a^{b}f(x)dx\simeq \sum _{i=1}^N{w}_{i}f(x_i)=\overrightarrow{w}\cdot \overrightarrow{f}}$$ とベクトル内積を使って評価してもよい。

（例題Pint－２）
下のプログラムは、ガウス・ルジャンドル積分によりN=4のガウス点で $${\int }_0^1{x}^{4}dx$$ の積分を計算するプログラムです。作成して実行しましょう。



gcc pexint-2.c -o pexint-2 -lm 
./pexint-2
      0.200000

この結果は解析値1/5と完全に一致している。これはN次のガウス・ルジャンドル積分が、N次多項式を用いて積分を評価しているため、f(x)=x^4のような多項式関数はN≧4で厳密値を与えてくれるためである。（N=1,2,3だと厳密値から外れる）
（課題Pint-３）
例題を参考にして、ガウス・ルジャンドル積分により f(x)=4/(1+x²) を a= 0 から b=1 まで積分するプログラムを作成し、実行しなさい。分割数を4, 8, 16, 32, 64, 128, 256, 512, ...としたとき、積分の結果の収束性を調べなさい。

（課題Pint-４)
ガウス・ルジャンドル積分により f(x)=exp(-x) についてa=0 かつ b=1の区間でガウス・ルジャンドル積分を実行するプログラムを作成しなさい。その結果と解析解とを比較せよ。

（余裕のある人へ）
（課題Pint-５)ガウス・ジャンドル積分と台形公式積分を比較し、どちらが優れているか？そして、なぜそちらの方が優れているか、その理由を考察しなさい。(実用上はガウス・ルジャンドル積分を使うことが多い)