概念の復習：広義積分

区間$[0,1]$における演習１の関数$f(x)=\frac{1}{\sqrt{x}}$ の積分を検討するとき、数値$f$は$x=0$で値が無限大になっていますので、 計算できないです。この場合、広義積分が考えられます。

【定義】以下の極限が存在するとき、この極限を$f$の$[a,b]$における広義積分と呼ぶ。

$$ \lim_{\epsilon \to 0+} \int_{a+\epsilon}^b f(x)~ dx $$

演習 2

以下の関数の広義積分を数値計算で検討してください。

$$ 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の値に対応する積分値のグラフを描いて、積分値が収束しているがとうかを観察します。

Step 1

関数f1とf2を定義します。また、積分値を計算する関数get_integralを作ります。今回の演習では、積分範囲$[a,1]$も変化しますので、 get_integralの変数にaを追加します。

In [21]:
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}) $$

In [22]:
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) );

f1の数値積分：1.366572 
正しい積分値（厳密値）: 1.367544 
f2の数値積分：2.298542 
正しい積分値（厳密値）: 2.302585

Step 2: aの値に応じて積分値を確認します。

まず、分割数をN=500として、f1の積分値を計算します。

In [23]:
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
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 1 1.5 2 2.5 3 3.5 4 4.5 5 gnuplot_plot_1a gnuplot_plot_2a

Step 3 観察

上記の計算結果によりますと、下積分が収束して、上積分が発散するそうです。数値積分の中の誤差の影響を確認するために、もっと大きい分割数で確認することが必要です。

N=1000で確認してみます。

In [24]:
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
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 1 1.5 2 2.5 3 3.5 4 4.5 5 gnuplot_plot_1a gnuplot_plot_2a

観察（続き）

上積分と下積分の差が小さくなっているので、N=10000で計算してみます。

In [30]:
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
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 1.2 1.4 1.6 1.8 2 2.2 1 1.5 2 2.5 3 3.5 4 4.5 5 gnuplot_plot_1a gnuplot_plot_2a

観察（続き）

N=10000のとき、上積分と下積分の差がもっと近づいています。 よって、分割数が十分大きいとき、上積分と下積分が一致することが見込まれます。

また、aの値が小さくなるほど、積分が一定になっていることが見られます。よって、広義積分が存在すると推定します。

Step 4

f2について、同じ方法で検討することが可能です。以下N＝50000のときの計算結果のみ説明します。

下のグラフによって、上積分と下積分両方が発散することを確認できます。よって、f2の$[0,1]$における広義積分が存在しないと推定します。

In [32]:
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
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 2 4 6 8 10 12 14 1 1.5 2 2.5 3 3.5 4 4.5 5 gnuplot_plot_1a gnuplot_plot_2a

結論

数値積分の計算結果によって、$[0,1]$では、

  • 関数$f_1$は広義積分を持って、広義積分の値が２になることを数値的に確認しました。
  • $f_2$が広義積分を持たないと考えられます。

補足

原始関数で計算した積分を利用して、以下のことも確認できます。

$$ \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; (発散する) $$