数値解析シラバス C言語環境と数値解析の概要(1 回) C言語環境と数値解析の概要(1 回) 方程式の根(1回) 方程式の根(1回) 連立一次方程式(2回) 連立一次方程式(2回) 補間と近似(2回) 補間と近似(2回) 数値積分(1回) 数値積分(1回) 中間試験(1回) 中間試験(1回) 常微分方程式(1回) 常微分方程式(1回) 偏微分方程式(2回) 偏微分方程式(2回) 数値計算の応用(3回) 数値計算の応用(3回) 期末試験(1回) 期末試験(1回) 6月12日(水)5限 今週 来週
数値積分 解析的に積分が困難な関数を数値的に積分する
数値積分 台形公式 台形公式 関数を直線で近似し、台形の集合体とす る。 関数を直線で近似し、台形の集合体とす る。 x y x0x0 x1x1 x2x2 x3x3 y0y0 y1y1 y2y2 y3y3
台形 I = 0.5×h {(y 0 +y 1 )+(y 1 +y 2 )+.. + (y n-1 +y n )} I = 0.5×h {(y 0 +y 1 )+(y 1 +y 2 )+.. + (y n-1 +y n )} = 0.5×h(y 0 +2y 1 + 2y y n ) = 0.5×h(y 0 +2y 1 + 2y y n ) x y x0x0 x1x1 x2x2 x3x3 y0y0 y1y1 y2y2 y3y3 h
台形公式計算例 x y I = x 4 +2x dx = I = x 4 +2x dx = I = 0.5×1×(0+2 ・ 3+2 ・ 20+87)=66.5 I = 0.5×1×(0+2 ・ 3+2 ・ 20+87)=66.5
Excelで確認 ABC 1XY 20.0=A2^4+2*A2=B2 30.5=2*B =B8 =(0.5/2)*sum(C2:C8) 区間の幅 数値の合計 I = x 4 +2x dx = I = x 4 +2x dx = N = 6
Excelで確認 ABC 1XY 20.0=A2^4+2*A2=B2 30.1=2*B =B8 =(0.1/2)*sum(C2:C32) 区間の幅 数値の合計 I = x 4 +2x dx = I = x 4 +2x dx = N = 30
台形公式による積分 #include #defineN30 doublefunc_y( double ); int main( ) { doubley[N+1]; doublexa = 0.0, xb = 3.0; doublez = 0.0, h, x, s; inti;
h = (xb - xa)/(double)N; for(i=0;i<=N;i++){ x = xa + h*(double)i; y[i] = func_y(x); } for(i=1;i<N;i++)z += 2.0*y[i]; s = (h/2.0)*(y[0]+z+y[N]); printf("ANS = %8.4lf\n",s); return 0; } doublefunc_y(double x){ returnpow(x, 4.0) + 2.0*x; } I =0.5×h(y 0 +2y 1 + 2y y n ) I =0.5×h(y 0 +2y 1 + 2y y n ) y = x 4 + 2x y = x 4 + 2x
数値積分での誤差縮小 分割数を多くする. 分割数を多くする. 曲線で近似 → シンプソン 曲線で近似 → シンプソン
厳密解の推定 nanswer ∞57.6
台形公式(第4章1)
台形公式(第4章1)
曲線による近似 シンプソン公式 シンプソン公式 2区間毎を二次関数で近似し、積分する。 x y x0x0 x1x1 x2x2 x3x3 y0y0 y1y1 y3y3 y4y4 y2y2 x4x4 2次関数
2次関数の積分 x y x0x0 x1x1 x2x2 y0y0 h y1y1 y2y2 h I = h(y 0 + 4y 1 + y 2 ) I = h(y 0 + 4y 1 + y 2 )13
2次関数の積分 x y y = x 2 I = x 2 dx =[ ] = I = x 2 dx =[ ] = 2 0 x3x333x3x 解析解 I = h(y 0 + 4y 1 + y 2 ) I = h(y 0 + 4y 1 + y 2 ) = 1 × (0 + 4 × 1 + 4)= = 1 × (0 + 4 × 1 + 4)= 数値解
シンプソン公式 x y x0x0 x1x1 x2x2 x3x3 y0y0 y1y1 y3y3 y4y4 x4x4 I = h{(y 0 + 4y 1 + y 2 )+ I = h{(y 0 + 4y 1 + y 2 )+ (y 2 + 4y 3 + y 4 ) + ……} (y 2 + 4y 3 + y 4 ) + ……}13 y2y2
シンプソン公式 x y x0x0 x1x1 x2x2 x3x3 y0y0 y1y1 y3y3 y4y4 x4x4 I = h{y 0 + 4(y 1 + y ) I = h{y 0 + 4(y 1 + y ) 2(y 2 + y ) + y n } 2(y 2 + y ) + y n }13 奇数番目を4 倍 偶数番目を2 倍 偶数番目 奇数番目 y2y2
Excelで確認 ABC 1XY 20.0=A2^4+2*A2=B2 30.5=?*B3 41.0=?*B4 51.5=?*B5 62.0=?*B6 72.5=?*B7 83.0=B8 =(?/?)*sum(C2:C8) 区間の幅数値の合計 I = x 4 +2x dx = I = x 4 +2x dx =
Excelで確認 ABC 1XY 20.0=A2^4+2*A2=B2 30.5=4*B3 41.0=2*B4 51.5=4*B5 62.0=2*B6 72.5=4*B7 83.0=B8 =(0.5/3)*sum(C2:C8) 初め 数値の合計 奇数番目 偶数番目 最後 I = x 4 +2x dx = I = x 4 +2x dx =
シンプソン計算例 x y I = x 4 +2x dx = I = x 4 +2x dx =
nclude #include #defineN30 doublefunc_y( double ); int main( ){ doubley[N+1]; doublexa = 0.0, xb = 3.0; doublez1 = 0.0, z2 = 0.0, h, x, s; inti; h = (xb - xa)/(double)N;
for(i=0;i<=N;i++){ x = xa + h*(double)i; y[i] = func_y(x); } for(i=1;i<=N-1;i+=2)z1 += 4.0*y[i]; for(i=2;i<=N-2;i+=2)z2 += 2.0*y[i]; s = (h/3.0)*(y[0]+z1+z2+y[N]); printf("ANS = %8.4lf\n",s); return 0; } doublefunc_y(double x) { returnpow(x, 4.0) + 2.0*x; } 奇数番目を4 倍 偶数番目を2 倍 最初と最後は1 倍
一次(台形)と二次(シンプソン) の精度 n台形 シンプソ ン ∞57.6
一次(台形)と二次(シンプソン) の精度