250 likes | 370 Views
数値解析シラバス. C言語環境と数値解析の概要(1回) 方程式の根(1回) 連立一次方程式(2回) 補間と近似(2回) 数値積分(1回) 中間試験(1回) 常微分方程式(1回) 偏微分方程式(2回) 数値計算の応用(3回) 期末試験(1回). 今週. 6月12日(水)5限. 来週. 数値積分. 解析的に積分が困難な関数を 数値的に積分する. 数値積分. 台形公式 関数を直線で近似し、台形の集合体とする。. y. y 3. y 1. y 2. y 0. x. x 0. x 1. x 2. x 3. 台形.
E N D
数値解析シラバス • C言語環境と数値解析の概要(1回) • 方程式の根(1回) • 連立一次方程式(2回) • 補間と近似(2回) • 数値積分(1回) • 中間試験(1回) • 常微分方程式(1回) • 偏微分方程式(2回) • 数値計算の応用(3回) • 期末試験(1回) 今週 6月12日(水)5限 来週
数値積分 解析的に積分が困難な関数を 数値的に積分する
数値積分 • 台形公式 関数を直線で近似し、台形の集合体とする。 y y3 y1 y2 y0 x x0 x1 x2 x3
台形 I = 0.5×h {(y0+y1)+(y1+y2)+ .. + (yn-1+yn)} = 0.5×h(y0+2y1 + 2y2 + .. + yn) y y3 y1 y2 y0 h x x0 x1 x2 x3
3 0 台形公式計算例 y I = 0.5×1×(0+2・3+2・20+87)=66.5 100 I = x4+2x dx= 57.6 75 50 25 x 3 0 1 2
3 0 Excelで確認 N = 6 I = x4+2x dx= 57.6 数値の合計 区間の幅
3 0 Excelで確認 N = 30 I = x4+2x dx= 57.6 数値の合計 区間の幅
台形公式による積分 #include <stdio.h> #include <math.h> #define N 30 double func_y( double ); int main( ) { double y[N+1]; double xa = 0.0, xb = 3.0; double z = 0.0, h, x, s; int i;
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; } double func_y(double x){ return pow(x, 4.0) + 2.0*x; } I =0.5×h(y0+2y1 + 2y2 + .. + yn) y = x4 + 2x
数値積分での誤差縮小 分割数を多くする. 曲線で近似 → シンプソン
台形公式(第4章1) 1 64 1 32 1 16 1 8 1 4
曲線による近似 • シンプソン公式 2区間毎を二次関数で近似し、積分する。 y 2次関数 2次関数 y4 y1 y3 y2 y0 x x0 x1 x2 x3 x4
1 3 I = h(y0 + 4y1 + y2) 2次関数の積分 y h h y2 y1 y0 x x0 x1 x2
数値解 4 8 3 x3 3 1 3 8 3 1 3 I = h(y0 + 4y1 + y2) = 1×(0 + 4×1 + 4)= y = x2 y 1 解析解 I = x2dx =[ ] = 2 0 0 1 2 x 2 0 2次関数の積分
1 3 シンプソン公式 I = h{(y0 + 4y1 + y2)+ (y2 + 4y3 + y4) + ……} y y4 y1 y3 y2 y0 x x0 x1 x2 x3 x4
1 3 奇数番目 シンプソン公式 奇数番目を4倍 I = h{y0 + 4(y1 + y3 + .. ) 2(y2 + y4 + .. ) + yn} 偶数番目を2倍 y 偶数番目 y4 y1 y3 y2 y0 x x0 x1 x2 x3 x4
3 0 Excelで確認 I = x4+2x dx= 57.6 区間の幅 数値の合計
3 0 Excelで確認 初め I = x4+2x dx= 57.6 偶数番目 奇数番目 数値の合計 最後
3 0 シンプソン計算例 y 32 I = x4+2x dx= 57.6 24 16 8 x 3 0 1 2
nclude <stdio.h> #include <math.h> #define N 30 double func_y( double ); int main( ){ double y[N+1]; double xa = 0.0, xb = 3.0; double z1 = 0.0, z2 = 0.0, h, x, s; int i; 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; } double func_y(double x) { return pow(x, 4.0) + 2.0*x; } 奇数番目を4倍 偶数番目を2倍 最初と最後は1倍
一次(台形)と二次(シンプソン)の精度 1 64 1 32 1 16 1 8 1 4