アルゴリズム演習2　 数値積分 岩崎　慶.

Presentation on theme: "アルゴリズム演習2　 数値積分 岩崎　慶."— Presentation transcript:

1 アルゴリズム演習2　 数値積分 岩崎　慶

2 演習内容 MATLABによる数値積分の習得 台形公式 シンプソン公式 integral, integral2, integral3

3 積分 区間[𝑎, 𝑏]で定義された関数𝑓(𝑥)の積分𝑆 𝑆= 𝑎 𝑏 𝑓 𝑥 𝑑𝑥 𝑓(𝑥) 𝑥 𝑎 𝑏 この面積と積分値Sは等しい

4 台形公式 小区間[ 𝑥 𝑖 , 𝑥 𝑖+1 ]の面積を台形の面積 𝑆 𝑖 で近似 𝑆= 𝑖=0 𝑛 𝑆 𝑖
𝑆= 𝑖=0 𝑛 𝑆 𝑖 𝑆 𝑖 = 1 2 (𝑓 𝑥 𝑖 +𝑓 𝑥 𝑖+1 )( 𝑥 𝑖+1 − 𝑥 𝑖 ) (𝑎= 𝑥 0 < 𝑥 1 <⋯< 𝑥 𝑛−1 < 𝑥 𝑛 =𝑏) 𝑥𝑖 𝑥 𝑖+1 𝑓( 𝑥 𝑖 ) 𝑓( 𝑥 𝑖+1 ) 𝑥 𝑆 𝑖 = 1 2 (𝑓 𝑥 𝑖 +𝑓 𝑥 𝑖+1 )( 𝑥 𝑖+1 − 𝑥 𝑖 )

5 台形公式 台形公式による積分の近似 等間隔に分割 𝑥 𝑖+1 − 𝑥 𝑖 =ℎとすると
𝑎 𝑏 𝑓 𝑥 𝑑𝑥≈ 𝑖=0 𝑛−1 𝑓 𝑥 𝑖 +𝑓 𝑥 𝑖 ( 𝑥 𝑖+1 − 𝑥 𝑖 ) (𝑎= 𝑥 0 < 𝑥 1 <⋯< 𝑥 𝑛−1 < 𝑥 𝑛 =𝑏) 𝑎 𝑏 𝑓 𝑥 𝑑𝑥≈ ℎ 2 𝑖=0 𝑛−1 (𝑓 𝑥 𝑖 +𝑓 𝑥 𝑖+1 ) = ℎ 2 (𝑓 𝑥 𝑓 𝑥 1 +⋯+𝑓 𝑥 𝑛−1 +𝑓 𝑥 𝑛 ) 台形公式

6 シンプソン公式 関数𝑓(𝑥)を2次多項式𝑝(𝑥)でラグランジュ近似
3点( 𝑥 0 ,𝑓( 𝑥 0 )),( 𝑥 1 ,𝑓( 𝑥 1 )),( 𝑥 2 ,𝑓( 𝑥 2 ))を𝑝(𝑥)が通る 3点が等間隔 𝑥 1 = 𝑥 0 +ℎ, 𝑥 2 = 𝑥 1 +ℎとすると 𝑝 𝑥 = 𝑎 0 𝑥− 𝑥 1 𝑥− 𝑥 2 + 𝑎 1 𝑥− 𝑥 0 𝑥− 𝑥 2 + 𝑎 2 𝑥− 𝑥 0 (𝑥− 𝑥 1 ) 𝑝 𝑥 0 = 𝑎 0 𝑥 0 − 𝑥 1 𝑥 0 − 𝑥 2 =𝑓 𝑥 0 𝑝 𝑥 1 = 𝑎 1 𝑥 1 − 𝑥 0 𝑥 1 − 𝑥 2 =𝑓 𝑥 1 𝑝 𝑥 2 = 𝑎 2 𝑥 2 − 𝑥 0 𝑥 2 − 𝑥 1 =𝑓( 𝑥 2 ) 𝑎 0 −ℎ −2ℎ =𝑓 𝑥 0 より 𝑎 0 = 𝑓 𝑥 ℎ 2 𝑎 1 ℎ −ℎ =𝑓 𝑥 1 より 𝑎 1 =− 𝑓 𝑥 1 ℎ 2 𝑎 2 2ℎ ℎ =𝑓 𝑥 2 より 𝑎 2 = 𝑓 𝑥 ℎ 2 𝑥0 𝑥 𝑓(𝑥0) 𝑓(𝑥1) 𝑓(𝑥2) 𝑥 1 𝑥 2 3点を通る2次多項式𝑝(𝑥)で𝑓(𝑥)を近似

7 シンプソン公式の導出 𝑥 0 𝑥 2 𝑝 𝑥 𝑑𝑥 = 𝑎 0 𝑥 0 𝑥 2 𝑥− 𝑥 1 𝑥− 𝑥 2 𝑑𝑥+ 𝑎 1 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 2 𝑑𝑥+ 𝑎 2 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 1 𝑑𝑥 𝐼 1 𝐼 2 𝐼 3 𝐼 1 = 𝑥 0 𝑥 2 𝑥− 𝑥 1 𝑥− 𝑥 2 𝑑𝑥= 𝑥 0 𝑥 2 𝑥− 𝑥 0 −ℎ 𝑥− 𝑥 0 −2ℎ 𝑑𝑥 𝑡=𝑥− 𝑥 0 −ℎ とおくと 𝑑𝑡=𝑑𝑥, 𝑥 : 𝑥 0 → 𝑥 0 +2ℎより𝑡:−ℎ→ℎ 𝐼 1 = 𝑥 0 𝑥 2 𝑥− 𝑥 0 −ℎ 𝑥− 𝑥 0 −2ℎ 𝑑𝑥 = −ℎ ℎ 𝑡 𝑡−ℎ 𝑑𝑡= −ℎ ℎ 𝑡 2 −ℎ𝑡 𝑑𝑡 = 𝑡 3 3 − ℎ𝑡 −ℎ ℎ = 2 ℎ 3 3

8 シンプソン公式の導出 𝑥 0 𝑥 2 𝑝 𝑥 𝑑𝑥 = 𝑎 0 𝑥 0 𝑥 2 𝑥− 𝑥 1 𝑥− 𝑥 2 𝑑𝑥+ 𝑎 1 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 2 𝑑𝑥+ 𝑎 2 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 1 𝑑𝑥 𝐼 1 𝐼 2 𝐼 3 𝐼 2 = 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 2 𝑑𝑥= 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 0 −2ℎ 𝑑𝑥 𝑡=𝑥− 𝑥 0 とおくと 𝑑𝑡=𝑑𝑥, 𝑥 : 𝑥 0 → 𝑥 0 +2ℎより𝑡:0→2ℎ 𝐼 2 = 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 0 −2ℎ 𝑑𝑥 = 0 2ℎ 𝑡 𝑡−2ℎ 𝑑𝑡= 0 2ℎ ( 𝑡 2 −2ℎ𝑡) 𝑑𝑡 = 𝑡 3 3 −ℎ 𝑡 ℎ =− 4 ℎ 3 3

9 シンプソン公式の導出 𝑥 0 𝑥 2 𝑝 𝑥 𝑑𝑥 = 𝑎 0 𝑥 0 𝑥 2 𝑥− 𝑥 1 𝑥− 𝑥 2 𝑑𝑥+ 𝑎 1 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 2 𝑑𝑥+ 𝑎 2 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 1 𝑑𝑥 𝐼 1 𝐼 2 𝐼 3 𝐼 3 = 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 1 𝑑𝑥= 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 0 −ℎ 𝑑𝑥 𝑡=𝑥− 𝑥 0 とおくと 𝑑𝑡=𝑑𝑥, 𝑥 : 𝑥 0 → 𝑥 0 +2ℎより𝑡:0→2ℎ 𝐼 3 = 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 0 −ℎ 𝑑𝑥 = 0 2ℎ 𝑡 𝑡−ℎ 𝑑𝑡= 0 2ℎ (𝑡 2 −ℎ𝑡) 𝑑𝑡 = 𝑡 3 3 − ℎ 𝑡 ℎ = 2 ℎ 3 3

10 シンプソン公式の導出 小区間[ 𝑥 0 𝑥 2 ]の𝑝(𝑥)の積分値
小区間[ 𝑥 0 𝑥 2 ]の𝑝(𝑥)の積分値 𝑥 0 𝑥 2 𝑝 𝑥 𝑑𝑥 = 𝑎 0 𝑥 0 𝑥 2 𝑥− 𝑥 1 𝑥− 𝑥 2 𝑑𝑥+ 𝑎 1 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 2 𝑑𝑥+ 𝑎 2 𝑥 0 𝑥 2 𝑥− 𝑥 0 𝑥− 𝑥 1 𝑑𝑥 𝐼 1 𝐼 2 𝐼 3 = 𝑎 0 2 ℎ 3 3 − 𝑎 1 4ℎ 𝑎 2 2 ℎ 3 3 = 𝑓 𝑥 ℎ ℎ 3 3 − − 𝑓 𝑥 1 ℎ ℎ 𝑓 𝑥 ℎ ℎ = ℎ 3 (𝑓 𝑥 0 +4𝑓 𝑥 1 +𝑓 𝑥 2 )

11 シンプソン公式 𝑎 = 𝑥0 <⋯< 𝑥𝑛 = 𝑏で𝑛を偶数とし、 𝑥 𝑖+1 − 𝑥 𝑖 =ℎ, (𝑖 = 0,…,𝑛−1)とすると 𝑥 0 𝑥 𝑛 𝑓 𝑥 𝑑𝑥≈ 𝑥 0 𝑥 2 𝑝 𝑥 𝑑𝑥 + 𝑥 2 𝑥 4 𝑝 𝑥 𝑑𝑥+⋯+ 𝑥 𝑛−2 𝑥 𝑛 𝑝 𝑥 𝑑𝑥 シンプソン公式

12 MATLABによる積分計算 1変数関数の積分 integral integral( f, xmin, xmax ) f : 被積分関数
f x ) sin( x ); I = integral( f, 0, pi / 2 ); >> I = integral( f, 0, pi/2 ) I = 1.0000

13 MATLABによる積分計算 2変数関数の積分 integral2
integral2( f, xmin, xmax, ymin, ymax ) f : 被積分関数 xmin, xmax, ymin, ymax : 積分区間 𝑥 𝑚𝑖𝑛 𝑥 𝑚𝑎𝑥 𝑦 𝑚𝑖𝑛 𝑦 𝑚𝑎𝑥 𝑓 𝑥,𝑦 𝑑𝑥𝑑𝑦 f x, y ) sin( x ) .* cos( y ); I0 = integral2( f, 0, pi / 2, -pi/2, pi/2 ); ymin x ) 1 – x; ymax x ) 1 + x; I1 = integral2( f, 0, pi/2, ymin, ymax ); >> I0 = integral2( f, 0, pi/2, -pi/2, pi/2 ) I0 = 2.0000 >> I1 = integral2( f, 0, pi/2, ymin, ymax ) I1 = 0.8487

14 レポートのヒント 台形公式による積分法の関数Mファイル % f : 被積分関数 % a , b : 積分区間 % n : 区間の分割数
function S = trap(f, a, b, n) S = f( a ) + f( b ); % S = f(x0)+f(xn) % forループで2*{f(x1)+f(x2)+....+f(xn-1)}をSに足す % h = (b – a) / n; 等間隔の分割h % xi = a + h * i; % Sにh/2=(b-a)/(2*n)を掛ける 台形公式の擬似プログラム

15 レポートのヒント シンプソン公式による積分法の関数Mファイル % f : 関数 % a , b : 積分区間
% n : 区間の分割数(nは偶数とする) function S = simpson(f, a, b, n) S = f( a ) + f( b ); %S = f(x0)+f(xn) % iが1からn-1まで % forループを偶数と奇数に分けて考える % iが偶数の場合 2*f(xi)をSに足す % iが奇数の場合 4*f(xi)をSに足す % Sにh/3=(b-a)/(3*n)を掛ける シンプソン公式の擬似プログラム