関数の積分法と収束性

区間$[a,b]$における関数$f$の積分を考えるとき、区間を分割して、分割点での関数の値を使って、関数の積分の近似値を求めることが可能です。 このような計算手法で得られる積分と数値積分を呼ばれています。

本日の授業では、区間の分割による計算した数値積分の収束性を検討します。即ち、分割の小区間のサイズが0に行くとき、積分値が発散するか、収束するか、ということを検討します。

1.数値積分の収束性

例として、前回と同じ関数の積分を考えます。

$$ \int_{0}^{1}{\sqrt {x}}+1dx = \frac{5}{3} $$

準備として、まず関数fを定義します。定義された関数は他のセルでしようできることを注意してください。

In [1]:
%関数fを定義する
function value=f(x)
    value = sqrt(x)+1;
end

以下のコードは与えられる分割のサイズによって、積分の上界と下界を計算しています。

$f$は単調増加しているので、各小区間$[x_1, x_2]$では、$f(x_1)$はこの小区間での$f$の下界、$f(x_2)$はこの小区間での$f$の上界を与えます。

In [60]:
h=0.1;  x=0:h:1;  n=length(x);
S1=0;
S2=0;
for k=1:(n-1)
  x1 = x(k);   y1 = f(x1); 
  x2 = x(k+1); y2 = f(x2);
  S1 = S1 + h*f(x1); %下界
  S2 = S2 + h*f(x2); %上界 
end
printf("h=%10.6f とき、 積分の下界= %15.10f , 積分の上界= %15.10f .\n", h, S1, S2)
h=  0.100000 とき、 積分の下界=    1.6105093417 , 積分の上界=    1.7105093417 .

上記のセルに書かれたコードを関数 get_integral(h)に定義します。

  • 入力: N (分割数)

  • 出力: [lower_bound, upper_bound] (積分値の下界、上界)

In [62]:
function [lower_bound, upper_bound] = get_integral(N)
    h = 1/ N;
    x=0:h:1;  n=length(x);
    S1=0;
    S2=0;
    for k=1:(n-1)
      x1 = x(k);   y1 = f(x1); 
      x2 = x(k+1); y2 = f(x2);
      S1 = S1 + h*f(x1); %下界
      S2 = S2 + h*f(x2); %上界 
    end
   lower_bound = S1;
   upper_bound = S2;
end

関数の戻り値をとるために、以下の形で関数を呼び出し、それぞれの戻り値がmylow, myupperに対応します。

In [45]:
[mylow, myupper] = get_integral(100)
mylow =  1.6615
myupper =  1.6715

積分値と分割のサイズhとの関係

分割のサイズ$h$がゼロに近づくとき、積分値の下界と上界の挙動を観察します。

実際の計算では、以下の分割数$N$と分割サイズ$h=1/N$の値を取って、下界と上界を計算して、それぞれを配列lower_list と upper_listに格納しなさい。 $$ N=4,8,16,32,64,128,256 $$

In [46]:
N_list = [4,8,16,32,64,128,256];
lower_list = [];
upper_list = [];

for N = N_list
  % 各hに対して、get_integralによって、積分の下界と上界を計算して、lower_listとuppper_listを更新します。
  %  ここにコードを書いてください。
  [mylow, myupper] = get_integral(N);
  lower_list = [lower_list, mylow];
  upper_list = [upper_list, myupper];  
end

lower_listとN_listとの関係図、upper_listとN_listとの関係図を描いてください。

In [47]:
hold on

plot(N_list, lower_list, 'ro-')
plot(N_list, upper_list, 'bo-')
grid on
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 1.5 1.55 1.6 1.65 1.7 1.75 1.8 0 50 100 150 200 250 300 gnuplot_plot_1a gnuplot_plot_2a

演習1:

  • 以下の関数について、$[0.1,1]$における数値積分の収束性を検討しなさい。

    $$ f(x)=\frac{1}{\sqrt{x}} $$

注意: get_integralの積分範囲を[0.1, 1]に修正して、分割数$N$に応じで、分割のサイズを再計算するのは必要です。

In [49]:
% 演習1に必要なコードをまとめて、ここに書いてください。

演習1の観察と結論:

(演習1の計算結果から分かることを書いてください。)

2. 広義積分

区間$[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 $$

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

In [59]:
a_list  =[0.1, 0.01, 0.001, 0.0001, 0.00001];
s_list = [];
for a = a_list
  % 各aに応じて、[a,1]における積分を計算しなさい。
  
  
end
plot(s_list, '-+')
grid on
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 "-"

演習2 の観察と結論:

(演習2の計算結果から分かることを書いてください。)

3. リーマン積分 (オプション)

区間$I=[a,b]$上の有界関数$f$について、分割$I$に対して、上積分と下積分の定義を考える。

まず、分割$I_h$を以下のように書く

$$ I=I_1 \cup I_2 \cup \cdots I_n $$

$h$は分割のサイズを示す。

$$ h = \max_i |I_i| \quad \quad (|I_i|: \mbox{区間}I_i\mbox{の幅}) $$

上積分と下積分は以下のように定義される。 $$ \mbox{上積分} =\sum_{i=1}^n I_i\mbox{における $f$ の最大値} \times |I_i| $$

$$ \mbox{下積分} =\sum_{i=1}^n I_i\mbox{における $f$ の最小値} \times |I_i| $$

リーマン積分の定義

【定義】以下の式が成り立つするとき、上積分と下積分の極限は$f$のリーマン積分と呼ぶ。($f$はリーマン積分の定義で可積分)

$$ \lim_{h\to 0} \mbox{上積分} = \lim_{h\to 0} \mbox{下積分} $$

演習 3

区間$I=[0,1]$とする。以下の関数はリーマン積分の定義で可積分ですか?適当な分割を選んで、上積分と下積分を計算しなさい。

特に、各分割の小区間での関数の単調性を確認する上、小区間における$f$の最小値と最大値を求めることが大切です。

  • $f(x)=sin(\pi x)$
  • $f(x)=sin(\pi x^2)$
In [43]:
x=0:0.001:1;
y1=sin(pi*x);
y2=sin(pi*x.^2);

plot(x,y1, x,y2)
grid on
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 gnuplot_plot_1a gnuplot_plot_2a
In [ ]: