連立一次方程式 a11x1+a12x2+a13x3+...+a1nxn= b1

Presentation on theme: "連立一次方程式 a11x1+a12x2+a13x3+...+a1nxn= b1"— Presentation transcript:

1 連立一次方程式 a11x1+a12x2+a13x3+...+a1nxn= b1 a21x1+a22x2+a23x3+...+a2nxn= b2 . an1x1+an2x2+an3x3+...+annxn= bn

2 身近な連立一次方程式 －ＣＴスキャン－

3 ＣＴスキャンの原理 未知数 ２５個 検査領域

4 ＣＴスキャンの原理 A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 A14 A15 A16 A17
X線管 X線センサ 検査領域

5 ＣＴスキャンの原理 A1+ A2 + A3 + A4 + A5 = 0 A6+ A7 + A8 + A9 + A10 = 3
検査領域

6 ＣＴスキャンの原理 検査領域 A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 A14 A15 A16
… … … A1 A2 A3 3 1 . A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 A14 A15 A16 A17 A18 A19 A20 A21 A22 A23 A24 A25 = （２５×２５）行列 … A25 A25 検査領域

7 実際のスキャンニング Ｘ線センサ Ｘ線管

8 身近な連立一次方程式 －ＣＴスキャン－ (3000×3000)×(3000×3000)行列 （９００万）×（９００万）行列 0.1 mmなら
解像度が 0.1 mmなら 少なくとも (3000×3000)×(3000×3000)行列 300 mm （９００万）×（９００万）行列 300 mm

9 連立方程式の解法の利用 シミュレーション 応力変形解析 流体解析 熱解析 有限要素法 境界要素法 差分法

10 有限要素法による構造解析

11 有限要素法による流体解析

12 掃き出し法 Sweep Out Method Gauss Elimination Gauss-Jordan
本日の解法 掃き出し法 Sweep Out Method Gauss Elimination Gauss-Jordan

13 連立一次方程式の解法 a00x0+a01x1+a02x2+...+a0n-1xn-1= b0 a10x1+a11x1+a12x2+...+a1n-1xn-1= b1 . an-10x0+an-11x2+...+an-1n-1xn-1= bn-1

14 教科書の例題 2x + y + 3z = 13 x + 3y + 2z = 13 3x + 2y + z = 10

15 紙とペンによる解法 2x + y + 3z = 13 x + 3y + 2z = 13 3x + 2y + z = 10 2 1 3 13

16 係数を１にする 2 1 3 13 10 消去する 1 0.5 1.5 6.5 2.5 -3.5 -9.5

17 係数を１にする 1 0.5 1.5 6.5 2.5 -3.5 -9.5 消去する 消去する 1 1.4 5.2 0.2 2.6 -3.6 -10.8

18 消去する 1 1.4 5.2 0.2 2.6 -3.6 -10.8 係数を１にする 1 2 3 答え

19 数式表示 a00 a01　 a a0n-1 x b0 a10 a11　 a a1n-1 x b1 = . . . an-10 an-11　an an-1n-1 xn bn-1

20 係数を１にする 式がどうなるか考えよう a00 a01　 a a0n b0 a10 a11　 a a1n b1 . . an-10 an-11　an an-1n bn-1

21 ピボットでわる a00/a00 a01/a00 a02/a00 …. a10 a11　 a a1n b1 消去する 式がどうなるか考えよう . . an-10 an-11　an an-1n bn-1

22 ピボットで割ったもの a01　 a a10-1*a a11-a01*a10 a12-a02*a10 ... 消去する . . この計算は黄色の行ではしない！！ an-10-1*an-10　an-11-a01*an-10...

23 消去する a a02 …. a a12 …. . . １にする 消去する an an

24 教科書 20P #include <stdio.h> #include <math.h> #define N 3
#define EPS .001 int main() { double a[N][N+1]={ { 2, 1, 3, 13}, { 1, 3, 2, 13}, { 3, 2, 1, 10} }; double pivot, del; int i, j, k, l; 教科書 20P

25 if( fabs(pivot) < EPS ){ printf("ピボット数が許容範囲以下\n"); return 1; }
for(i=0;i<N;i++){ pivot = a[i][i]; if( fabs(pivot) < EPS ){ printf("ピボット数が許容範囲以下\n"); return 1; } for( j=i; j < N+1; j++) a[i][j]=a[i][j] / pivot; for(k=0; k<N;k++){ del = a[k][i]; for( j=i;j<N+1; j++) if( k!=i) a[k][j]-=del*a[i][j]; 対角行を代入 ゼロ割を避ける 対角行を１にする 等しくないとき真 代入演算子

26 2 1 3 13 10 for(i=0;i<N;i++){ pivot = a[i][i];
if( fabs(pivot) < EPS ){ printf("ピボット数が許容範囲以下\n"); return 1; } for( j=i; j < N+1; j++) a[i][j]=a[i][j] / pivot; for(k=0; k<N;k++){ del = a[k][i]; for( j=i;j<N+1; j++) if( k!=i) a[k][j]-=del*a[i][j]; 対角行を代入 １にする １にする 2 1 3 13 10

27 1 0.5 1.5 6.5 3 2 13 10 for(i=0;i<N;i++){ pivot = a[i][i];
if( fabs(pivot) < EPS ){ printf("ピボット数が許容範囲以下\n"); return 1; } for( j=i; j < N+1; j++) a[i][j]=a[i][j] / pivot; for(k=0; k<N;k++){ del = a[k][i]; for( j=i;j<N+1; j++) if( k!=i) a[k][j]-=del*a[i][j]; 1 0.5 1.5 6.5 3 2 13 10 ０にする １にする

28 printf("X%d = %6.2lf\n", l, a[l][N] ); return 0; }
for( l=0; l<N; l++) printf("X%d = %6.2lf\n", l, a[l][N] ); return 0; } 出力

29 教科書の結果

30 演習問題2.1 3x + y + z = 10 x + 5y + 2z = 21 x + 2y + 5z = 30

31 演習2.1 #include <stdio.h> #include <math.h> #define N 3
#define EPS .001 int main() { double a[N][N+1]={ { 3, 1, 1, 10}, { 1, 5, 2, 21}, { 1, 2, 5, 30} }; double pivot, del; int i, j, k, l; 演習2.1 変更

32 演習問題2.1の結果