4. 連立一次方程式の解法 4. LU 分解法 同じ係数行列 A( サイズ n n) をもつ m 組の連立 次方程式 AX = B ( ただし A=[ ij ] は n 行 n 列の正則行列,B=[b ij ] と X=[x ij ] は n 行 m 列の行列 ) を同時に解く 行列 A,B を並置して 個の配列 A (n 行 n+m 列 ) を作成し, i,n+j =b ij (i=,,n; j=,,m) として A の配列成分を改めて ij と記す A' = [ A B] = M n n L O L n n 以下のアルゴリズムを実行すると, M nn b b b M n b b b n L O L b m b m = M M bnm n n A' の B 部に解 X の値が入る L O L n n M nn, n+, n+ M n, n+, n+, n+ n, m+ L O L, n+ m, n+ m M n, n+ m アルゴリズム 4.. LU 分解法の素直なアルゴリズム LU 分解と前進消去の同時進行 ( 原型 ): for j= to n+m {k= {L: l i = i j := j / {U: u j = j /l {Y: y j =b j /l for k= to n {k=, 3,, n for i=k to n {L: l ik = ik Σ k- l= l il u lk for l= to k- ik := ik - il lk for j=k+ to n+m {U: u kj = ( kj -Σ k- l= l kl u lj )/l kk for l= to k- {Y: y kj = (b kj -Σ k- l= l kl y lj )/l kk kj := kj - kl lj kj := kj / kk 後退代入 : for k=n-, n-,, {k=n {X: x nj = y nj for j= to m { 右辺ベクトルの数だけ繰り返す for l=k+ to n {k n {X: x kj =y kj Σ n l=k+ u kl x lj k,n+j := k,n+j - kl l,n+j プログラム 4.. 3 4 /* Solve Liner Eqution System AX=B */ /* By LU Decomposition Method */ /* ( Originl Algorithm ) */ /* AX = B */
5 6 7 8 9 0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 36 37 38 39 40 4 4 43 44 45 46 47 48 49 50 5 5 53 54 55 /* A : Coefficient Mtrix size n*n */ /* X : Solution Mtrix (m Vectors) size n*m */ /* B : Right-hnd-side Mtrix (m Vectors) size n*m */ /* */ /* HERE A'=[A,B] */ # include <stdio.h> #define MMAX 0 #define NMAX 30 void min(void) { flot [NMAX][NMAX+MMAX]; int n,m; int i,j,k,l; /* Red Dt */ printf("input Mtrix Size: n, m? "); /* 行列データの画面入力 */ scnf("%d%d", &n, &m); /* サイズ :n, m */ printf(" ninput Coefficient Mtrix A n"); /* 行列 A の画面入力 */ printf("? A(%d, j), j=,%d : ", i+, n); for(j=0; j<n; j++) scnf("%g", &[i][j]); printf(" ninput R.H.S. Vectors B n"); /* 行列 B の画面入力 */ printf("? B(%d, j), j=,%d : ", i+, m); for(j=0; j<m; j++) scnf("%g", &[i][n+j]); /* Write Dt */ printf(" nliner Eqution System: AX=B n"); printf("* Size of Mtrix: n,m=%d,%d n", n,m); printf("* A'=[A B] n"); /* 行列 A'=[A B] の画面への出力 */ for(j=0;j<n+m;j++) printf("%6.g", [i][j]); /* LU Decomposition nd Forwrd reduction */ for(j=;j<n+m;j++) [0][j] /= [0][0]; for(k=;k<n;k++){ for(i=k;i<n;i++) for(l=0;l<k;l++) [i][k] -= [i][l]*[l][k]; for(j=k+;j<n+m;j++){
56 57 58 59 60 6 6 63 64 65 66 67 68 69 70 7 7 73 74 75 76 77 78 79 80 8 8 83 84 85 for(l=0;l<k;l++) [k][j] -= [k][l]*[l][j]; [k][j] /= [k][k]; printf(" nmtrix A' fter LU decomposition nd Forwrd reduction n"); /* Debug Write */ for(j=0;j<n+m;j++) printf("%5.6e", [i][j]); /* Bckwrd Substitution */ for(k=n-;k>=0;--k){ for(j=0;j<m;j++){ for(l=k+;l<n;l++){ [k][n+j] -= [k][l]*[l][n+j]; /* Print Result */ printf(" nsolutions n"); for(j=0;j<m;j++) printf("%5.6e", [i][n+j]); アルゴリズム 4.. LU 分解法のコンパクトなアルゴリズム LU 分解と前進消去の同時進行 : for k= to n for j=k+ to n+m kj := kj / kk {U( ただし対角成分 を除く ) for i=k+ to n ij := ij - ik kj {L, U と Y 後退代入 : for k=n, n-,, for j= to m for i= to k- i,n+j := i,n+j - ik k,n+j 3
次のプログラムでは, 入力データをファイルから読んでいる プログラム 4.. 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 36 37 38 39 40 4 4 43 44 45 /* Solve Liner Eqution System AX=B */ /* By LU Decomposition Method */ /* ( Compct Algorithm ) */ /* AX = B */ /* A : Coefficient Mtrix size n*n */ /* X : Solution Mtrix (m Vectors) size n*m */ /* B : Right-hnd-side Mtrix (m Vectors) size n*m */ /* */ /* HERE A'=[A,B] */ /* */ /* Input Dt is red from File */ # include <stdio.h> #define MMAX 0 #define NMAX 30 flot [NMAX][NMAX+MMAX]; int n,m; void input(void) /* 関数 :A, B の係数をファイルから入力 */ { int i,j; FILE *fp; /* ファイルポインターの宣言 */ /* Red Dt */ /* ファイルオープンとエラー処理 */ if( (fp=fopen("c4mtrix.dt","r"))==null ){ printf(" ファイルをオープンできません n"); exit(); /* 行列データのファイル入力 */ fscnf(fp,"%d%d",&n, &m); /* サイズ :n, m */ /* 行列 A */ for(j=0; j<n; j++) fscnf(fp,"%g", &[i][j]); /* 行列 B */ for(j=0; j<m; j++) fscnf(fp,"%g", &[i][n+j]); /* ファイルクローズ */ fclose(fp); /* 行列データの出力 */ printf("liner Eqution System: AX=B n"); 4
46 47 48 49 50 5 5 53 54 55 56 57 58 59 60 6 6 63 64 65 66 67 68 69 70 7 7 73 74 75 76 77 78 79 80 8 8 83 84 85 86 87 88 89 90 9 9 93 94 printf("* Size of Mtrix: n,m=%d,%d n", n,m); printf("* A'=[A B] n"); /* 行列 A'=[A B] の画面への出力 */ for(j=0;j<n+m;j++) printf("%6.g", [i][j]); void min(void) /* メイン関数 */ { int i,j,k,l; input(); /* LU Decomposition nd Forwrd reduction */ for(k=0;k<n;k++){ for(j=k+;j<n+m;j++){ [k][j] /= [k][k]; for(i=k+;i<n;i++){ [i][j] -= [i][k]*[k][j]; printf(" nmtrix A' fter LU decomposition nd Forwrd reduction n"); /* Debug Write */ for(j=0;j<n+m;j++) printf("%5.6e", [i][j]); /* Bckwrd Substitution */ for(k=n-;k>=0;--k){ for(j=0;j<m;j++){ for(i=0;i<k;i++){ [i][n+j] -= [i][k]*[k][n+j]; /* Print Result */ printf(" nsolutions n"); for(j=0;j<m;j++) printf("%5.6e", [i][n+j]); 5
実行例 問題 ) デバッグ ( プログラムの誤りを修正すること ) 用 問題 ) Ax = b 3 6 A = 7 b = 4 0 5 AX = B 5 3 0 0 0 5 3 0 0 0 A = B = 3 5 0 0 0 4 5 8 0 0 0 問題 に対する実行例 ) プログラム 4..( 入力データをキーボードから入力 ) ------------------------- 実行開始 ----------------------- Input Mtrix Size: n, m? 3 Input Coefficient Mtrix A? A(, j), j=, 3 : 3? A(, j), j=, 3 : 7? A( 3, j), j=, 3 : 4 0 Input R.H.S. Vectors B? B(, j), j=, : 6? B(, j), j=, :? B( 3, j), j=, : 5 Liner Eqution System: AX=B * Size of Mtrix: n,m=3, * A'=[A B] 3 6 7 4 0 5 Mtrix A' fter LU decomposition nd Forwrd reduction.000000e+00.000000e+00 3.000000e+00 6.000000e+00.000000e+00 3.000000e+00.000000e+00 3.000000e+00.000000e+00.000000e+00 3.000000e+00.000000e+00 Solutions.000000e+00.000000e+00.000000e+00 ------------------------- おしまい ----------------------- 6
問題 に対する実行例 ) プログラム 4..( 入力データをファイルから読む ) 入力データ ( ファイル名 :c4mtrix.dt): ------------------------- 入力データ ----------------------- 4 5 5 3-5 - 3-4 5 0 0 0 3 0 0 0 5 0 0 0-8 0 0 0 -------------------------------------------------------- ------------------------- 実行開始 ----------------------- Liner Eqution System: AX=B * Size of Mtrix: n,m=4,5 * A'=[A B] 5 3-0 0 0 5 3 0 0 0-3 - 5 0 0 0 4 5-8 0 0 0 Mtrix A' fter LU decomposition nd Forwrd reduction.000000e+00.500000e+00.500000e+00-5.000000e-0.000000e+00 5.000000e-0 0.000000e+00 0.000000e+00 0.000000e+00 5.000000e+00 -.050000e+0 6.90476e-0-3.333333e-0.90476e-0.38095e-0-9.5380e-0 0.000000e+00 0.000000e+00.000000e+00-4.500000e+00 4.8574e+00-4.666667e-0 3.466667e+00.333333e-0 -.000000e-0.333333e-0 0.000000e+00.000000e+00 -.000000e+00 -.38095e+00 5.0e+00 -.000000e+00 -.50443e-0-4.64608e-0 6.45930e-0.995e-0 Solutions.000000e+00-5.309733e-0.477876e-0-8.84956e-03-6.9469e-0 -.000000e+00.50445e-0-3.53983e-0 -.4599e-0 8.849559e-03 3.000000e+00 7.96460e-0 -.684e-0.63743e-0 9.9036e-0 -.000000e+00 -.50443e-0-4.64608e-0 6.45930e-0.995e-0 ------------------------- おしまい ----------------------- 7
4. ガウスの消去法 4. 節と同様, B AX = を解くにあたり, 配列を A '= [ A B] にとる 4.. ガウスの消去法 ( 基本 ) アルゴリズム 4.. 前進消去 : for k= to n- for i=k+ to n w:= ik/ kk for j=k+ to n+m {j=k は ik =0 と自明なので省く ij := ij - w kj { (k+) ij= (k) ij ( (k) ik (k) kk) (k) kj {b (k+) ij= b (k) ij ( (k) ik (k) kk)b (k) kj 後退代入 : for k=n, n-,, for j= to m { 右辺ベクトルの数だけ繰り返す if k<n then {k<n: x kj =(b (k) kj Σ n l=k+ (k) kl x lj )/ (k) kk for l=k+ to n {k=n: x kj = b (k) kj/ (k) kk k,n+j := k,n+j - kl l,n+j k,n+j := k,n+j / kk 後退代入はつぎのようにコンパクトなアルゴリズムに書き換えられる : for k=n, n-,, for j= to m { 右辺ベクトルの数だけ繰り返す k,n+j := k,n+j/ k,k {x kj =(b (k) kj Σ n l=k+ (k) kl x lj )/ (k) kk for i= to k- i,n+j := i,n+j - ik k,n+j プログラム 4.. LU 分解法のプログラム 4.. に記した関数 input を用い, ファイルからデータを入力する 以下にガウスの消去法 ( 後退代入にはコンパクトなアルゴリズムを使用 ) の min 関数を記す 他はプログラム 4.. と同じ 3 4 5 void min(void) /* メイン関数 */ { int i,j,k; flot w; 8
6 7 8 9 0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 input(); /* Forwrd reduction */ for(k=0;k<n-;k++){ for(i=k+;i<n;i++){ w=[i][k]/[k][k]; for(j=k+;j<n+m;j++){ [i][j] -= w*[k][j]; /* Bckwrd Substitution */ for(k=n-;k>=0;--k){ for(j=0;j<m;j++){ [k][n+j] /= [k][k]; for(i=0;i<=k-;i++){ [i][n+j] -= [i][k]*[k][n+j]; /* Print Result */ printf(" nsolutions n"); for(j=0;j<m;j++) printf("%5.6e", [i][n+j]); 実行例 教科書の例題 4. を解く 4 4 Ax = b A = 3 4 b = 7 3 8 8 入力データ ------------------------- 入力データ ----------------------- 3 4 3 4 3 8 4 7 8 -------------------------------------------------------- ------------------------- 実行開始 ----------------------- Liner Eqution System: AX=B * Size of Mtrix: n,m=3, * A'=[A B] 4 4 9
3 4 7 3 8 8 Solutions -3.000000e+00.000000e+00.000000e+00 ------------------------- おしまい ----------------------- 4.. 掃出し法 ( あるいはガウス ジョルダン消去法 ) アルゴリズム 4.. 掃き出しのプロセス : for k= to n w:=/ kk for j=k to n+m { (k+) kj= (k) kj/ (k) kk kj := kj w {b (k+) kj= b (k) kj/ (k) kk for i= to n if(i k) then w:= ik for j=k to n+m ij := ij - w kj { (k+) ij= (k) ij (k) ik (k+) kj {b (k+) ij= b (k) ij (k) ik b (k+) kj プログラム 4.. LU 分解法のプログラム 4.. に記した関数 input を用い, ファイルからデータを入力する 以下に min 関数を記す ( 他はプログラム 4.. と同じ ) 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 0 void min(void) /* メイン関数 */ { int i,j,k; flot w; input(); /* Sweep-out */ for(k=0;k<n;k++){ w = /[k][k]; for(j=k;j<(n+m);j++) [k][j] *= w; if(i!=k){ w = [i][k]; for(j=k;j<(n+m);j++) [i][j] -= w*[k][j]; 0
3 4 5 6 7 8 9 /* Print Result */ printf(" nsolutions n"); for(j=0;j<m;j++) printf("%5.6e", [i][n+j]); 実行例 LU 分解法の実行例で記した問題 を解く ------------------------- 実行開始 ----------------------- Liner Eqution System: AX=B * Size of Mtrix: n,m=4,5 * A'=[A B] 5 3-0 0 0 5 3 0 0 0-3 - 5 0 0 0 4 5-8 0 0 0 Solutions.000000e+00-5.309733e-0.477876e-0-8.849546e-03-6.9469e-0 -.000000e+00.50445e-0-3.53983e-0 -.4599e-0 8.849557e-03 3.000000e+00 7.96460e-0 -.684e-0.63743e-0 9.9036e-0 -.000000e+00 -.50443e-0-4.64608e-0 6.45930e-0.995e-0 -------------------------おしまい----------------------- 4..3 ピボット選択付ガウス消去法 ( あるいはピボット選択付 LU 分解法 ( ただし L の対角成分はすべて )) アルゴリズム 4..3 LU 分解と前進消去 : for i= to n p i =i { 置換情報 p の初期設定 for k= to n- i=k,k+, n に対し, ik > kk となる r を探す if r k then {k 行と r 行の p,l,a,b の値を入れ換える p k と p r の値を入れ換える k 行の kj と r 行の rj の値を入れ換える (j=,,,n+m) for i=k+ to n ik := ik/ kk {L( ただし対角成分 を除く ) for j=k+ to n+m ij := ij - ik kj {L,U と Y 後退代入は, アルゴリズム 4.. におけるものと同じ
プログラム 4..3 LU 分解法のプログラム 4.. に記した関数 input を用い, ファイルからデータを入力する 以下に min 関数を記す ( 他はプログラム 4.. と同じ ) 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 36 37 38 39 40 4 4 43 44 45 46 47 48 49 50 5 5 void min(void) /* メイン関数 */ { int p[nmax],i,j,k,r; flot mx,w,kk,ptemp,temp; input(); for(i=0;i<n;i++) /* Index for Pivoting */ p[i]=i; /* Forwrd Substitution with Pivoting */ for(k=0;k<n-;k++){ mx=fbs([k][k]); r=k; for(i=k+;i<n;i++){ /* serch r such tht _{r,k=mx _{i,k i=k,..n */ if( fbs([i][k]) > mx ){ mx=fbs([i][k]); r=i; if(r!=k){ /* exchnge column r nd column k */ ptemp=p[r]; p[r] =p[k]; p[k] =ptemp; for(j=0;j<n+m;j++){ temp = [r][j]; [r][j] = [k][j]; [k][j] = temp; for(i=k+;i<n;i++) [i][k] = [i][k]/[k][k]; for(i=k+;i<n;i++){ for(j=k+;j<n+m;j++) [i][j] -= [i][k]*[k][j]; /* Bckwrd Substitution */ for(k=n-;k>=0;--k){ kk=/[k][k]; for(j=0;j<m;j++){ for(i=k+;i<n;i++){ [k][n+j] -= [k][i]*[i][n+j]; [k][n+j] *=kk; /* print result */ printf(" npermuttion n"); printf("%3d", p[i]); if(i%0==9)
53 54 55 56 57 58 59 60 6 6 63 64 65 66 printf(" nlu decomposition n"); for(j=0;j<n;j++) printf("%5.6e", [i][j]); printf(" nsolution n"); for(j=0;j<m;j++) printf("%5.6e", [i][n+j]); 実行例 ガウスの消去法 ( 基本 ) における実行例で扱った問題 ( 教科書例題 4.) を解く ------------------------- 実行開始 ----------------------- Liner Eqution System: AX=B * Size of Mtrix: n,m=3, * A'=[A B] 4 4 3 4 7 3 8 8 Permuttion 0 LU decomposition 3.000000e+00 8.000000e+00.00000e+0 6.666667e-0 -.333333e+00-5.333333e+00 3.333333e-0 -.499999e-0-9.999996e-0 Solution -3.000000e+00.000000e+00 9.999998e-0 ------------------------- おしまい ----------------------- 3
4.3 三項方程式の解法 三項方程式 c i x i- + i x i + b i x i+ = d i (i=,,,n) は, つぎのように行列表示される Ax = d b x d c b x d O O O M M A =, x =, d = O O O M M c n n b n x n dn cn n xn d これは 4.,4. 節に示した数値解法により解くことができる つぎのアルゴリズムは 4. 節の LU 分解によるものであり, これを実行すると右辺ベクトルの d 部に解 x の値が入る なお,c, b n の値は使用しないが, 通常 c =0, b n =0 を入れておく アルゴリズム 4.3. LU 分解と前進消去の同時進行 : b :=b / d :=d / for i= to n i := i -c i b i- d i :=(d i -c i h i- )/ i b i :=b i / i 後退代入 : for i=n-, n-,, d i :=d i -b i d i+ 問題 常微分方程式 d u/dx = を条件 u(0)=0, u()=0 のもとで解くと, 解 u(x)=0.5x(x-) を得るが, これを差分法で求め, 差分近似解と誤差を出力する x の区間 [0,] を n- 等分して各小区間幅を Δx とし,x i =(i-)δx,i=,,,n,( ただし Δx=/(n-)) とおく 0 x x x 3 x n- x n Δx 階導関数は u"(x i ) ( u(x i +Δx) u(x i ) + u(x i -Δx) )/(Δx) と近似される ( 第 9 章 : 数値微分を参照のこと ) ので,U i =u(x i ) と記せば,i=,3,,n- に対し U i+ U i + U i- = (Δx) が成立し,i= と i=n に対しては以下の条件が課せられる U = 0 U n = 0 n 個の未知数 U i, i=,,,n に対し, 代数方程式が n 個得られたので, 連立一次方程式により U i が求められる 具体的に代数式を列挙すると 4
x=x =0 のとき U = 0 x=x のとき U - U + U 3 = (Δx) x=x 3 のとき U U 3 + U 4 = (Δx) x=x 4 のとき U 3 U 4 + U 5 = (Δx) : : : : x=x n- のとき U n- U n- + U n = (Δx) x=x n = のとき U n = 0 となり, 三項方程式の形をしていることがわかる プログラム 4.3. 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 36 37 38 39 40 4 4 43 /* three-term equtions */ #include <stdio.h> #define NMAX 0 void min(void) { flot [NMAX],b[NMAX],c[NMAX],d[NMAX], xmin=0.0,xmx=.0,dx,x,error; int n,i; printf("number of points: n? "); scnf("%d", &n); dx=(xmx-xmin)/(n-); c[0]= 0.0; [0]=.0; b[0]= 0.0; d[0]= 0.0; for(i=;i<n-;i++){ c[i]=.0; [i]=-.0; b[i]=.0; d[i]= dx*dx; c[n-]= 0.0; [n-]=.0; b[n-]= 0.0; d[n-]= 0.0; 5 /* dt */ printf(" nthree-term equtions: n i,c[i],[i],b[i],d[i] n"); /* output */ for(i=0;i<n;i++) printf("%d %g %g %g %g n", i,c[i], [i], b[i], d[i]); /* forwrd */ b[0] /= [0]; d[0] /= [0]; for(i=;i<n;i++){ [i] -= c[i]*b[i-]; d[i] -= c[i]*d[i-]; d[i] /= [i]; b[i] /= [i];
44 45 46 47 48 49 50 5 5 53 54 for(i=n-;i>=0;i--) d[i] -=b[i]*d[i+]; printf(" n i tx t Solution t Error n"); x=i*dx; error=0.5*x*(x-)-d[i]; printf("%d t%5.6g t%5.6g t%5.6e n", i,x,d[i],error); /* bckwrd */ 実行例 ------------------------- 実行開始 ----------------------- 近似解のグラフ number of points: n? Three-Term equtions: i,c[i],[i],b[i],d[i] 0 0 0 0-0.005-0.005 3-0.005 4-0.005 5-0.005 : : 8-0.005 9-0.005 0 0 0 0 Solution 0 0. 0.4 0.6 0.8 0-0.05-0. -0.5 i x Solution Error 0 0 0 0.000000e+00 0.05-0.0375 -.495944e-09 0. -0.045-6.58488e-09 3 0.5-0.06375-9.44984e-09 4 0. -0.08 -.758337e-08 5 0.5-0.09375 -.38306e-08 6 0.3-0.05 -.65407e-08 7 0.35-0.375 -.985580e-08 8 0.4-0. -.8794e-08 9 0.45-0.375 -.4949e-08 0 0.5-0.5 -.3574e-08 0.55-0.375 -.074987e-08 0.6-0. -.668930e-08 3 0.65-0.375 -.07004e-08 4 0.7-0.05-8.64673e-09 5 0.75-0.09375-4.65663e-09 6 0.8-0.08-5.6644e-09 7 0.85-0.06375-4.09578e-09 8 0.9-0.045-4.0333e-09 9 0.95\ -0.0375.346933e-09 0 0 7.45058e-09 -------------------------おしまい----------------------- 6
4.4 反復法 4.4. ヤコビ法 アルゴリズム 4.4. for i= to n { 初期値の設定 x i := 与えられた初期値 for k=0,,, { 収束するまで反復する for i= to n x temp :=b i for j= to n if(i j) then x temp :=x temp - ij x i x temp :=x temp / ii x new i:=x temp ===== for i= to n x i :=x new i ===== 収束判定は ノルムを用いて行い, 解の変化量が次式を満たすとき収束したとする n ( k+ ) ( k) ( k+ ) ( k) xi xi x i= ( k + ) n x ( k + ) xi i= x = < ε アルゴリズムの のところに分子 分母のノルムの計算 ( 下記 ) を挿入する norm_dx:=0 norm_x :=0 for i= to n norm_dx:=norm_dx + x new i-xi norm_x :=norm_x + x new i あるいは上記を, 直前の (for i= to n) のループと同時進行させる アルゴリズムの のところに次の収束判定を入れる If norm_dx/norm_x < ε brek(for k=0,,, のループから出る ) プログラム 4.4. 以下に示すプログラムでは, 関数 input においてファイルからデータを入力する 3 4 5 /* Solve Liner Eqution System Ax=b */ /* by Jcobi Method */ #include <stdio.h> #include <mth.h> 7
6 7 8 9 0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 36 37 38 39 40 4 4 43 44 45 46 47 48 49 50 5 5 53 54 55 56 57 58 59 60 6 6 #define KMAX 50 #define NMAX 0 flot eps, [NMAX][NMAX], b[nmax]; int n; void input(void) /* 入力データを読み込む関数 */ { int i, j; FILE *fp; /* ファイルポインターの宣言 */ /* ファイルオープンとエラー処理 */ if( (fp=fopen("c4iter.dt","r"))==null ){ printf(" ファイルをオープンできません n"); exit(); fscnf(fp,"%g %d", &eps, &n); /* 収束判定許容誤差と行列サイズの入力 */ /* 行列データのファイル入力 */ for(j=0;j<n;j++) fscnf(fp,"%g",&[i][j]); for(i=0;i<n;i++) /* 右辺ベクトルのファイル入力 */ fscnf(fp,"%g",&b[i]); fclose(fp); /* ファイルクローズ */ printf("solve Liner Eqution System: Ax=b n"); printf("* Mtrix size: n=%d t Tolernce: eps=%g n n",n,eps); /* 画面への出力 */ printf("* A'=[A b] n"); /* 行列 A'=[A b] の画面への出力 */ for(j=0;j<n;j++) printf("%5.6e", [i][j]); printf("%5.6e n", b[i]); void min(void) /* メイン関数 */ { int i, j, k; flot dxnorm, xnorm,xx,x[nmax],xnew[nmax]; input(); printf("itertion by Jcobi Method n"); /* initil vlue */ for(i=0;i<n;i++) x[i]=0; /* itertion */ for(k=0;k<kmax;k++){ printf("k=%d t",k); for(i=0;i<n;i++) printf("%5.6e", x[i]); dxnorm=0; xnorm =0; 8
63 64 65 66 67 68 69 70 7 7 73 74 75 76 77 78 79 80 8 8 xx=b[i]; for(j=0;j<n;j++) if(i!=j) xx -= [i][j]*x[j]; xx /= [i][i]; dxnorm += fbs(xx-x[i]); xnorm += fbs(xx); xnew[i]= xx; for(i=0;i<n;i++) x[i]=xnew[i]; if(dxnorm/xnorm<eps) brek; if(k>=kmax) printf("not convergent? n"); /* output solution */ printf(" nsolution:"); for(i=0;i<n;i++) printf("%5.6g",x[i]); 4.4. ガウス ザイデル法 アルゴリズム 4.4. for i= to n { 初期値の設定 x i := 与えられた初期値 for k=0,,, { 収束するまで反復する for i= to n x temp :=b i for j= to n if(i j) then x temp :=x temp - ij x i x temp :=x temp / ii x i :=x temp プログラム 4.4. ヤコビ法のプログラム 4.4. に記した関数 input を用い, ファイルからデータを入力する 以下に min 関数を記す ( 他はプログラム 4.4. と同じ ) 3 4 5 6 7 8 9 void min(void) /* メイン関数 */ { int i, j, k; flot dxnorm, xnorm,xx,x[nmax]; input(); printf("itertion by Guss-Seidel Method n"); /* initil vlue */ 9
0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 36 37 38 for(i=0;i<n;i++) x[i]=0; /* itertion */ for(k=0;k<kmax;k++){ printf("k=%d t",k); for(i=0;i<n;i++) printf("%5.6e", x[i]); dxnorm=0; xnorm =0; xx=b[i]; for(j=0;j<n;j++) if(i!=j) xx -= [i][j]*x[j]; xx /= [i][i]; dxnorm += fbs(xx-x[i]); xnorm += fbs(xx); x[i]=xx; if(dxnorm/xnorm<eps) brek; if(k>=kmax) printf("not convergent? n"); /* output solution */ printf(" nsolution:"); for(i=0;i<n;i++) printf("%5.6g",x[i]); 4.4.3 SOR 法 ( 加速緩和法 ) アルゴリズム 4.4.3 for i= to n { 初期値の設定 x i := 与えられた初期値 for k=0,,, { 収束するまで反復する for i= to n x temp :=b i for j= to n if(i j) then x temp :=x temp - ij x i x temp :=x temp / ii x i :=x i + ω(x temp - x i ) プログラム 4.4.3 ヤコビ法のプログラム 4.4. に記した関数 input を用い, ファイルからデータを入力する 以下に min 関数を記す ( 他はプログラム 4.4. と同じ ) 0
3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 36 37 38 39 40 4 void min(void) /* メイン関数 */ { int i, j, k; flot dxnorm, xnorm,xx,dx,x[nmax]; flot omeg=.; input(); printf("itertion by SOR n"); /* initil vlue */ for(i=0;i<n;i++) x[i]=0; /* itertion */ for(k=0;k<kmax;k++){ printf("k=%d t",k); for(i=0;i<n;i++) printf("%5.6e", x[i]); dxnorm=0; xnorm =0; xx=b[i]; for(j=0;j<n;j++) if(i!=j) xx-= [i][j]*x[j]; xx /= [i][i]; dx = xx-x[i]; xx = x[i] + omeg*dx; dxnorm += fbs(dx); xnorm += fbs(xx); x[i]=xx; if(dxnorm/xnorm<eps) brek; if(k>=kmax) printf("not convergent? n"); /* output solution */ printf(" nsolution:"); for(i=0;i<n;i++) printf("%5.6g",x[i]); 実行例 次の問題を, ヤコビ法, ガウス ザイデル法,SOR 法で解け 5 9 5 6 Ax = b A = b = 5 4 4 入力データ ( ファイル名 :c4iter.dt): ------------------------- 入力データ -----------------------.0e-6
4-5 - - 5 5 - -4 9 6 - -4 -------------------------------------------------------- 実行例 ) ヤコビ法 ------------------------- 実行開始 ----------------------- Solve Liner Eqution System: Ax=b * Mtrix size: n=4 Tolernce: eps=e-06 * A'=[A b] -5.000000e+00.000000e+00 -.000000e+00.000000e+00 9.000000e+00 -.000000e+00 5.000000e+00.000000e+00.000000e+00 6.000000e+00.000000e+00.000000e+00 5.000000e+00.000000e+00 -.000000e+00.000000e+00.000000e+00 -.000000e+00-4.000000e+00-4.000000e+00 Itertion by Jcobi Method k=0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 k= -.800000e+00.00000e+00-4.000000e-0.000000e+00 k= -.00000e+00 7.00000e-0-9.00000e-0 5.00000e-0 k=3 -.88000e+00.044000e+00-6.480000e-0 8.00000e-0 k=4 -.70000e+00 9.300000e-0-9.040000e-0 8.90000e-0 k=5 -.086600e+00 9.80000e-0-8.696000e-0 8.735000e-0 k=6 -.0860e+00 9.89000e-0-9.44800e-0 9.93500e-0 k=7 -.049958e+00 9.847740e-0-9.44480e-0 9.359650e-0 k=8 -.03853e+00 9.96650e-0-9.583040e-0 9.57765e-0 k=9 -.06890e+00 9.95749e-0-9.79460e-0 9.68458e-0 k=0 -.0903e+00 9.95380e-0-9.79083e-0 9.77685e-0 : : k=34 -.000006e+00 9.999985e-0-9.999936e-0 9.99993e-0 k=35 -.000004e+00 9.999989e-0-9.999955e-0 9.99995e-0 k=36 -.000003e+00 9.99999e-0-9.999967e-0 9.999964e-0 Solution: - 0.999999-0.999998 0.999997 ------------------------- おしまい ----------------------- 実行例 ) ガウス ザイデル法 ------------------------- 実行開始 ----------------------- Solve Liner Eqution System: Ax=b * Mtrix size: n=4 Tolernce: eps=e-06 * A'=[A b] -5.000000e+00.000000e+00 -.000000e+00.000000e+00 9.000000e+00 -.000000e+00 5.000000e+00.000000e+00.000000e+00 6.000000e+00.000000e+00.000000e+00 5.000000e+00.000000e+00 -.000000e+00.000000e+00.000000e+00 -.000000e+00-4.000000e+00-4.000000e+00 Itertion by Guss-Seidel Method
k=0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 k= -.800000e+00 8.400000e-0-3.760000e-0 4.040000e-0 k= -.400800e+00 9.4400e-0-6.47360e-0 6.899440e-0 k=3 -.0309e+00 9.473767e-0-8.08665e-0 8.94064e-0 k=4 -.097e+00 9.737e-0-8.965963e-0 9.070589e-0 k=5 -.065503e+00 9.848068e-0-9.436457e-0 9.49366e-0 k=6 -.035708e+00 9.975e-0-9.6989e-0 9.73970e-0 k=7 -.0946e+00 9.954860e-0-9.83609e-0 9.849558e-0 k=8 -.00607e+00 9.975396e-0-9.908767e-0 9.98004e-0 k=9 -.00578e+00 9.986590e-0-9.95075e-0 9.95530e-0 k=0 -.0035e+00 9.9969e-0-9.97898e-0 9.975643e-0 : : k=0 -.000007e+00 9.999983e-0-9.999937e-0 9.999944e-0 k= -.000004e+00 9.99999e-0-9.999966e-0 9.999970e-0 k= -.00000e+00 9.999995e-0-9.99998e-0 9.999983e-0 Solution: - -0.999999 0.999999 ------------------------- おしまい ----------------------- 実行例 3)SOR 法 : 収束加速パラメータ ω=. の場合 ------------------------- 実行開始 ----------------------- Solve Liner Eqution System: Ax=b * Mtrix size: n=4 Tolernce: eps=e-06 * A'=[A b] -5.000000e+00.000000e+00 -.000000e+00.000000e+00 9.000000e+00 -.000000e+00 5.000000e+00.000000e+00.000000e+00 6.000000e+00.000000e+00.000000e+00 5.000000e+00.000000e+00 -.000000e+00.000000e+00.000000e+00 -.000000e+00-4.000000e+00-4.000000e+00 Itertion by SOR k=0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 k= -.60000e+00 9.6000e-0-4.039680e-0 3.06704e-0 k= -.4050e+00 9.85090e-0-7.78900e-0 9.054794e-0 k=3 -.494e+00 9.30703e-0-9.505594e-0 9.48769e-0 k=4 -.037866e+00.036e+00-9.66307e-0 9.88760e-0 k=5 -.00854e+00 9.9080e-0-9.943387e-0 9.9749e-0 k=6 -.00578e+00.005e+00-9.969443e-0 9.977609e-0 k=7 -.00070e+00 9.9948e-0-9.990863e-0 9.995770e-0 k=8 -.00054e+00 9.998699e-0-9.997874e-0 9.996573e-0 k=9 -.00007e+00.00003e+00-9.998674e-0 9.999738e-0 k=0 -.00004e+00 9.999583e-0-9.99984e-0 9.999635e-0 k= -.00008e+00.000009e+00-9.999856e-0 9.999947e-0 k= -.00000e+00 9.999955e-0-9.999976e-0 9.999976e-0 k=3 -.00000e+00.000000e+00-9.999989e-0 9.999989e-0 Solution: - - ------------------------- おしまい ----------------------- 参考文献 葉子 : 数値計算の基礎 解法と誤差, コロナ社 (007) 森口繁一, 伊理正夫, 武市正人編 :C による算法痛論, 東京大学出版会 (000) Heth, Michel T.: Scientific Computing, An Introductory Survey, McGrw-Hill(00) 3