区間$[0,1]$における演習1の関数$f(x)=\frac{1}{\sqrt{x}}$ の積分を検討するとき、数値$f$は$x=0$で値が無限大になっていますので、 計算できないです。この場合、広義積分が考えられます。
【定義】以下の極限が存在するとき、この極限を$f$の$[a,b]$における広義積分と呼ぶ。
$$ \lim_{\epsilon \to 0+} \int_{a+\epsilon}^b f(x)~ dx $$
以下の関数の広義積分を数値計算で検討してください。
$$ f_1(x)=\frac{1}{\sqrt{x}}, \quad f_2(x)=\frac{1}{x} $$
数値積分を計算を行うとき、分割数を$N=500$に固定して、積分範囲$[a,1]$の左辺$a$を以下の値を取ってください。 $$ a = 0.1, 0.01, 0.001, 0.0001, 0.00001 $$
即ち、各aの値に対応する積分値のグラフを描いて、積分値が収束しているがとうかを観察します。
関数f1とf2を定義します。また、積分値を計算する関数get_integralを作ります。今回の演習では、積分範囲$[a,1]$も変化しますので、 get_integralの変数にaを追加します。
function y=f1(x)
y= 1/sqrt(x);
end
function y=f2(x)
y= 1/x;
end
function [lower_bound, upper_bound] = get_integral_f1(a,N)
h = (1-a)/ N;
x=a:h:1; n=length(x);
S1=0;
S2=0;
for k=1:(n-1)
x1 = x(k);
x2 = x(k+1);
S1 = S1 + h*f1(x1); % 単調減少関数なので、S1が上界となる。
S2 = S2 + h*f1(x2); % 単調減少関数なので、S1が下界となる。
end
lower_bound = S2;
upper_bound = S1;
end
function [lower_bound, upper_bound] = get_integral_f2(a,N)
h = (1-a)/ N;
x=a:h:1; n=length(x);
S1=0;
S2=0;
for k=1:(n-1)
x1 = x(k);
x2 = x(k+1);
S1 = S1 + h*f2(x1); % 単調減少関数なので、S1が上界となる。
S2 = S2 + h*f2(x2); % 単調減少関数なので、S1が下界となる。
end
lower_bound = S2;
upper_bound = S1;
end
念のために、定義した関数を使って積分を計算してみます。
関数f1とf2の原始関数を求めて、積分関数の厳密値を計算します。数値積分で計算した値と厳密値を比較します。
$$ \int_a^b f_1(x) dx = 2 (\sqrt{b}- \sqrt{a}), \quad \int_a^b f_2(x) dx = (\log{b}- \log{a}) $$
v1 = get_integral_f1(0.1,1000);
printf("f1の数値積分:%f \n", v1);
printf("正しい積分値(厳密値): %f \n", 2*(sqrt(1) - sqrt(0.1)) );
v2=get_integral_f2(0.1,1000);
printf("f2の数値積分:%f \n", v2);
printf("正しい積分値(厳密値): %f \n", log(1) - log(0.1) );
まず、分割数をN=500として、f1の積分値を計算します。
a_list =[0.1, 0.01, 0.001, 0.0001, 0.00001];
s_list_low = [];
s_list_up = [];
for a = a_list
[v1,v2] = get_integral_f1(a, 500);
s_list_low = [s_list_low, v1];
s_list_up = [s_list_up, v2];
end
hold on
plot(s_list_low, '-+')
plot(s_list_up, '-+')
grid on
上記の計算結果によりますと、下積分が収束して、上積分が発散するそうです。数値積分の中の誤差の影響を確認するために、もっと大きい分割数で確認することが必要です。
N=1000で確認してみます。
a_list =[0.1, 0.01, 0.001, 0.0001, 0.00001];
s_list_low = [];
s_list_up = [];
for a = a_list
[v1,v2] = get_integral_f1(a, 1000);
s_list_low = [s_list_low, v1];
s_list_up = [s_list_up, v2];
end
hold on
plot(s_list_low, '-+')
plot(s_list_up, '-+')
grid on
上積分と下積分の差が小さくなっているので、N=10000で計算してみます。
a_list =[0.1, 0.01, 0.001, 0.0001, 0.00001];
s_list_low = [];
s_list_up = [];
for a = a_list
[v1,v2] = get_integral_f1(a, 10000);
s_list_low = [s_list_low, v1];
s_list_up = [s_list_up, v2];
end
hold on
plot(s_list_low, 'b-+')
plot(s_list_up, 'r-+')
grid on
N=10000のとき、上積分と下積分の差がもっと近づいています。 よって、分割数が十分大きいとき、上積分と下積分が一致することが見込まれます。
また、aの値が小さくなるほど、積分が一定になっていることが見られます。よって、広義積分が存在すると推定します。
f2について、同じ方法で検討することが可能です。以下N=50000のときの計算結果のみ説明します。
下のグラフによって、上積分と下積分両方が発散することを確認できます。よって、f2の$[0,1]$における広義積分が存在しないと推定します。
a_list =[0.1, 0.01, 0.001, 0.0001, 0.00001];
s_list_low = [];
s_list_up = [];
for a = a_list
[v1,v2] = get_integral_f2(a, 50000);
s_list_low = [s_list_low, v1];
s_list_up = [s_list_up, v2];
end
hold on
plot(s_list_low, 'b-+')
plot(s_list_up, 'r-+')
grid on
原始関数で計算した積分を利用して、以下のことも確認できます。
$$ \lim_{a\to 0+} \int_{a}^1 f_1(x) dx = \lim_{a\to 0+} 2(1-\sqrt{a}) = 2; (収束いる) $$ $$ \lim_{a\to 0+} \int_{a}^1 f_2(x) dx = \lim_{a\to 0+} (-\log{a}) = \infty; (発散する) $$