精度保証数値計算演習3

学籍番号： 600p037-6 長友泰崇

演習 3. 数値積分命令を使って、次の定積分の値を計算せよ。

(1) ∫_{0}^{1} e^{x} cos (x)dx  (2) ∫0^{\frac{π}{2}} \frac{1}{\sqrt{1 − 0.3 sin^{2} x}}dx
(3) ∫_{0}^{2π} \frac{sin (x)}{x}dx  (4) ∫_{0}^{1} \frac{e^{−x}}{\sqrt{x}}dx
(5) ∫_{0}^{0.2} \frac{1 − e^{−x}}{\sqrt{x}} dx
(6) (1)〜(5) の積分の 1 つを選び、その誤差を厳密に評価せよ。

1 環境

実行環境は以下の通りである。

OS : Vine Linux 2.0 CPU : Mobile Intel Celeron 500MHz

数値計算ソフト : Scilab-2.5

2 (1)〜(5)

(1)〜(5) を Scilab の数値積分命令によって以下のように積分の値を求める。

(1)
-->integrate('exp(x)*cos(x)','x',0,1)
 ans  =
    1.3780246
(2)
-->integrate('1/sqrt(1-0.3*(sin(x)^2))','x',0,%pi/2)
 ans  =
    1.7138894
(3)
-->integrate('sin(x)/x','x',0,2*%pi)

 ans  =
    1.4181516
(4)
-->integrate('exp(-x)/sqrt(x)','x',0,1)
 ans  =
    1.4936483
(5)
-->integrate('(1-exp(-x))/sqrt(x)','x',0,0.2)
 ans  =
     .0562147

3 (6)

3.1 積分の誤差評価

上で求めた

(1) ∫_{0}^{1} e^{x} cos (x)dx

について厳密な誤差評価を行う。

まず、積分の区間包囲を求める。関数 f を連続な関数とする定積分

I = ∫_{a}^{b} f(t)dt

を考える。区間 [a, b] 上に分点

t_{0} = a < t_{1} < t_{2} < ⋯ < t_{n} = b

をとる。 Δ_{i} = [t_{i−1}, t_{i}] とする。区間 [a, b] で考える単関数 u = ∑ u_{i} を分割 Δ = {Δ_{1}, Δ_{2}, ⋯, Δ_{n}} 上の単関数ということにする。

u(t) ≤ f(t) をみたす分割 Δ = {Δ_{1}, Δ_{2}, ⋯, Δ_{n}} 上の単関数 u(t) と、 v(t) ≥ f(t) をみたす分割 Δ′ = {Δ_{1}, Δ_{2}, ⋯, Δ_{n}} 上の単関数 v(t) があれば

∫_{a}^{b} u(t)dt ≤ I ≤ ∫_{a}^{b} v(t)dt

が成立する。

上の式の最左辺の積分の値の上限値である下ダルブー積分 I_{l} と最右辺の積分の値の下限値である上ダルブー積分 I_{u} の値が一致した値を I とするのであるが、上の式の最右辺および最左辺は短冊の面積の総和として簡単に求めることができる。したがって、次のプログラムは I の値を精度保証付きで求めるプログラムとなる。

call(′down′);

I_{l} = ∑_{i=1}^{n} u_{i}d(Δ_{i});

call(′up′);

I_{u} = ∑_{i=1}^{n} v_{i}d(Δ′_{i});

ただし、 d(Δ_{i}) と d(Δ′_{i}) はそれぞれ区間 Δ_{i} と Δ′_{i} の幅とする。このプログラムによって計算された値を用いると、

I ∈ [I_{l}, I_{u}]

となることがわかる。

次に誤差評価である。

t_{i} − t_{i−1} = h, h = (b − a)/n とする。

F_{i} ∈ f([t_{i−1}, t_{i}])

であれば、次の誤差評価を与える

∥ I − ∑_{k=1}^{n} F_{k}h ∥_{C[a,b]} ≤ h max_{t∈[a,b]} |f′(t)|

が成立する。

ここで、

F_{i} ∈ f([t_{i−1}, t_{i}]) = [min(f(t_{i−1}), f(t_{i})) , max(f(t_{i−1}), f(t_{i}))]

より、

∑_{k=1}^{n} F_{k}h ∈ [I_{l}, I_{u}]

となる。これより、

∥ I − ∑_{k=1}^{n} F_{k}h ∥_{C[a,b]} ⊂ ∥ I − [I_{l}, I_{u}] ∥_{C[a,b]}

であるから、

∥ I − [I_{l}, I_{u}] ∥_{C[a,b]} ≤ h max_{t∈[a,b]} |f′(t)|

が成り立つ。ここで、 Scilab の積分命令で求めた値を I_{s} とすると、

∥ I_{s} − [I_{l}, I_{u}] ∥_{C[a,b]} ≤ h max_{t∈[a,b]} |f′(t)|

が成り立てば、 Scilab で求めた値 I_{s} の誤差評価は妥当である、といえる。

4 プログラム

” scilab.sci ”に、以下のように記述する。

function  [I_l,I_u,I_s] = sekibun(a,b,n)
h=(b-a)/n;
call('down');
i = a;
I_l=0;
for j=a:n-1 ,
I_l=I_l+min(exp(i)*cos(i),exp(i+h)*cos(i+h))*h ,     //min(f(t_{i-1}),f(t_i))
i=i+h , end;
call('up');
i = a;
I_u=0;
for j=a:n-1 ,
I_u=I_u+max(exp(i)*cos(i),exp(i+h)*cos(i+h))*h ,     //max(f(t_{i-1}),f(t_i))
i=i+h ,  end;
printf(" I= [%lg , %lg]",I_l,I_u);                    //I=[I_l,I_u]
I_s=integrate('exp(x)*cos(x)','x',0,1);    //Scilab での数値積分における値:I_s
k = 1;                                         //k=max|f'(t)|
x=[abs(I_s-I_l) ; abs(I_s-I_u)];
if norm(x,"inf") <= h*k then,
     printf("%lg <= %lg ",norm(x,"inf"),h*k)   //norm(x,"inf")<=h*k
     printf("Errors evalution is O.K.!!"),     //誤差評価→ O.K. !!
else,

       printf("Errors evaluation is not good!!"), end
endfuntion

この関数は、 a, b, n を呼び込んで、 I ∈ [I_{l}, I_{u}] となる I_{u} と I_{u}、 (1) を Scilab の積分命令により求めた I_{s} を表示し、 I_{s} の誤差評価を行い、その評価が妥当ならば、

∥ I_{s} − [I_{l}, I_{u}] ∥_{C[a,b]} ≤ h max_{t∈[a,b]} |f′(t)|

を表示した上で、

”Errors evalution is O.K.!!”

と表示する関数である。このプログラム中にでてくる k は、

k = max_{t∈[a,b]} |f′(t)|

で、 (1) では、

k = max_{t∈[a,b]} |f′(t)| = max_{t∈[0,1]} |{e^{x} cos (x)}′| = max_{t∈[0,1]} |e^{x}{cos (x) − sin (x)}| = 1

となっている。

5 実行結果

前述した” scilab.sci ”を用いて、実際に (1) を Scilab の積分命令で求めた値 I_{s} の誤差評価を、 n = 10000 として行うと以下のようになる。

-->link('down.o','down','C');
linking files down.o  to create a shared executable
shared archive loaded
Linking down (in fact down)
Link done
-->link('up.o','up','C');
linking files up.o  to create a shared executable
shared archive loaded
Linking up (in fact up)
Link done
-->getf('sekibun.sci');
-->a=0
 a  =
    0.
-->b=1
 b  =

-->n=10000
n  =
10000.
-->[I_l,I_u,I_s]=sekibun(a,b,n)
I= [1.37799 , 1.37806]
3.16551e-05 <= 0.0001
Errors evalution is O.K.!!
I_s  =
1.3780246
I_u  =
1.3780563
I_l  =
1.377993

この結果より、 (1) について、 Scilab の積分命令で求められた値 I_{s} の誤差評価は妥当であることが証明された。

参考文献

[1] 大石進一 : ”応用解析セミナー数値計算”, 裳華房,1999.