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

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

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

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

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

ＣＴスキャンの原理 検査領域 A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 A14 A15 A16 1 1 1 1 1 0 0 0 0 0 … 0 0 0 0 0 1 1 1 1 1 0 … 0 0 0 0 0 0 0 0 0 0 1 1 1 … 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 = （２５×２５）行列 … 0 0 0 0 0 1 1 1 1 1 A25 A25 検査領域

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

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

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

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

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

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

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

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

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

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

係数を１にする 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

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

数式表示 a00 a01　 a02 . . . a0n-1 x0 b0 a10 a11　 a12 . . . a1n-1 x1 b1 = . . . an-10 an-11　an-12 . . an-1n-1 xn-1 bn-1

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

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

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

消去する 1 a01 a02 …. 0 a11 a12 …. . . １にする 消去する 0 an-11 an-12 ....

教科書 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

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]; 対角行を代入 ゼロ割を避ける 対角行を１にする 等しくないとき真 代入演算子

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

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 ０にする １にする

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; } 出力

教科書の結果

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

演習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 変更

演習問題2.1の結果