C言語による数値計算プログラミング演習

Size: px

Start display at page:

Download "C言語による数値計算プログラミング演習"

えのすみい
5 years ago
Views:

1 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 プログラム /* Solve Liner Eqution System AX=B */ /* By LU Decomposition Method */ /* ( Originl Algorithm ) */ /* AX = B */

2 /* 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++){

3 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 次のプログラムでは, 入力データをファイルから読んでいるプログラム /* 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

5 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

6 実行例問題 ) デバッグ ( プログラムの誤りを修正すること ) 用問題 ) Ax = b 3 6 A = 7 b = AX = B A = B = 問題に対する実行例 ) プログラム 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] Mtrix A' fter LU decomposition nd Forwrd reduction e e e e e e e e e e e e+00 Solutions e e e おしまい

7 問題に対する実行例 ) プログラム 4..( 入力データをファイルから読む ) 入力データ ( ファイル名 :c4mtrix.dt): 入力データ実行開始 Liner Eqution System: AX=B * Size of Mtrix: n,m=4,5 * A'=[A B] Mtrix A' fter LU decomposition nd Forwrd reduction e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e-0.995e-0 Solutions e e e e e e e e e e e e e e e e e e e-0.995e おしまい

8 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.. と同じ void min(void) /* メイン関数 */ { int i,j,k; flot w; 8

9 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 = 入力データ入力データ実行開始 Liner Eqution System: AX=B * Size of Mtrix: n,m=3, * A'=[A B] 4 4 9

10 Solutions e e e おしまい掃出し法 ( あるいはガウスジョルダン消去法 ) アルゴリズム 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.. と同じ ) 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

11 /* 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] Solutions e e e e e e e e e e e e e e e e e e e-0.995e おしまいピボット選択付ガウス消去法 ( あるいはピボット選択付 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.. におけるものと同じ

12 プログラム 4..3 LU 分解法のプログラム 4.. に記した関数 input を用い, ファイルからデータを入力する以下に min 関数を記す ( 他はプログラム 4.. と同じ ) 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)

13 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] Permuttion 0 LU decomposition e e e e e e e e e-0 Solution e e e おしまい

14 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

15 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 となり, 三項方程式の形をしていることがわかるプログラム /* 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];

16 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] : : Solution i x Solution Error e e e e e e e e e e e e e e e e e e e \ e e おしまい

17 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 においてファイルからデータを入力する /* Solve Liner Eqution System Ax=b */ /* by Jcobi Method */ #include <stdio.h> #include <mth.h> 7

18 #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

19 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. と同じ ) void min(void) /* メイン関数 */ { int i, j, k; flot dxnorm, xnorm,xx,x[nmax]; input(); printf("itertion by Guss-Seidel Method n"); /* initil vlue */ 9

20 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]); SOR 法 ( 加速緩和法 ) アルゴリズム 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. に記した関数 input を用い, ファイルからデータを入力する以下に min 関数を記す ( 他はプログラム 4.4. と同じ ) 0

21 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 法で解け Ax = b A = b = 入力データ ( ファイル名 :c4iter.dt): 入力データ e-6

22 実行例 ) ヤコビ法実行開始 Solve Liner Eqution System: Ax=b * Mtrix size: n=4 Tolernce: eps=e-06 * A'=[A b] e e e e e e e e e e e e e e e e e e e e+00 Itertion by Jcobi Method k= e e e e+00 k= e e e e+00 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 : : k= e e e e-0 k= e e e e-0 k= e e e e-0 Solution: おしまい実行例 ) ガウスザイデル法実行開始 Solve Liner Eqution System: Ax=b * Mtrix size: n=4 Tolernce: eps=e-06 * A'=[A b] e e e e e e e e e e e e e e e e e e e e+00 Itertion by Guss-Seidel Method

23 k= e e e e+00 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 : : k= e e e e-0 k= e e e e-0 k= e e e e-0 Solution: おしまい実行例 3)SOR 法 : 収束加速パラメータ ω=. の場合実行開始 Solve Liner Eqution System: Ax=b * Mtrix size: n=4 Tolernce: eps=e-06 * A'=[A b] e e e e e e e e e e e e e e e e e e e e+00 Itertion by SOR k= e e e e+00 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 k= e e e e-0 Solution: おしまい参考文献葉子 : 数値計算の基礎解法と誤差, コロナ社 (007) 森口繁一, 伊理正夫, 武市正人編 :C による算法痛論, 東京大学出版会 (000) Heth, Michel T.: Scientific Computing, An Introductory Survey, McGrw-Hill(00) 3

Ｃ言語による数値計算プログラミング演習

Ｃ言語による数値計算プログラミング演習 8. 数値積分任意の区間 [,b] における f() の定積分 b () I = f ( ) d の値は, つぎのように n 点の関数値の和により近似的に与えられる () In = Ak f ( k) n k = このとき, k を分点,A k を重みという 8. ニュートンコーツ ( 複合型 ) 積分公式積分区間 [,b] を等分割して n 個の分点をとり, 被積分関数 f() を n- 次ラグランジュ補間多項式で近似して得られる積分公式を