320 likes | 467 Views
連立一次方程式. a 11 x 1 +a 12 x 2 +a 13 x 3 +...+a 1n x n = b 1. a 21 x 1 +a 22 x 2 +a 23 x 3 +...+a 2n x n = b 2. a n1 x 1 +a n2 x 2 +a n3 x 3 +...+a nn x n = b n. 身近な連立一次方程式 -CTスキャン-. CTスキャンの原理. 未知数 25個. 検査領域. X 線管. X 線センサ. CTスキャンの原理. 検査領域. CTスキャンの原理. A 1 + A 2 + A 3 + A 4 + A 5 = 0.
E N D
連立一次方程式 a11x1+a12x2+a13x3+...+a1nxn= b1 a21x1+a22x2+a23x3+...+a2nxn= b2 . . . an1x1+an2x2+an3x3+...+annxn= bn
CTスキャンの原理 未知数 25個 検査領域
X線管 X線センサ CTスキャンの原理 検査領域
CTスキャンの原理 A1+ A2+ A3+ A4+ A5= 0 A6+ A7+ A8+ A9+ A10= 3 A11+ A12+ A13+ A14+ A15= 1 A16+ A17+ A18+ A19+ A20= 1 A21+ A22+ A23+ A24+ A25= 0 検査領域
CTスキャンの原理 A1 A2 A3 0 3 1 . . . 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 … = (25×25)行列 … 0 0 0 0 0 1 1 1 1 1 A25 A25 検査領域
実際のスキャンニング X線センサ X線管
身近な連立一次方程式-CTスキャン- 解像度が 0.1 mmなら 少なくとも (3000×3000)×(3000×3000)行列 300 mm (900万)×(900万)行列 300 mm
連立方程式の解法の利用 • 有限要素法 • 境界要素法 • 差分法 シミュレーション • 応力変形解析 • 流体解析 • 熱解析
本日の解法 掃き出し法 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
係数を1にする 消去する
係数を1にする 消去する 消去する
消去する 係数を1にする 答え
数式表示 a00 a01a02 . . . a0n-1x0 b0 a10 a11a12 . . . a1n-1x1 b1 = . . . . . . . . . an-10 an-11an-12. . an-1n-1xn-1 bn-1
係数を1にする 式がどうなるか考えよう a00 a01a02 . . . a0n-1 b0 a10 a11a12 . . . a1n-1 b1 . . . . . . an-10 an-11an-12. . an-1n-1 bn-1
ピボットでわる a00/a00a01/a00a02/a00 …. a10 a11a12 . . . a1n-1 b1 消去する . . . . . . an-10 an-11an-12. . an-1n-1 bn-1 式がどうなるか考えよう
ピボットで割ったもの 1 a01a02 . . . a10-1*a10 a11-a01*a10a12-a02*a10 ... 消去する . . . . . . この計算は黄色の行ではしない!! an-10-1*an-10an-11-a01*an-10...
消去する 1 a01a02 …. 0 a11 a12 …. . . . . . . 1にする 消去する 0 an-11 an-12 ....
#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
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にする 等しくないとき真 代入演算子
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にする 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]; } } 1にする 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
#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 変更