２０２３年前期 数値解析 第 9 回講義ノート
２０２３年６月１５日 三浦憲二郎
スライド１ 講義アウトライン•

連立方程式の反復解法•「•

反復法の原理：「•「•「」行列と「」行列」の原理」法」法

スライド２ 連立 1 次方程式 p.93•

行列のほとんどの要素が 0「•」行列

逆にそのほとんどが 0 でない「•sparse matrix」行列

dense matrix

疎行列を効率よく解く方法「」解法

スライド３ 反復法の原理
連立 1 次方程式
𝑨𝑨𝒙𝒙=𝒃𝒃(5.1)

を同値な方程式
𝒙𝒙=𝝓𝝓(𝒙𝒙)=𝑴𝑴𝑴𝑴+𝑵𝑵𝑵𝑵に変形し，「」法により解を求める．𝑴𝑴を「「」行列という．

」写像の原理スライド 4 縮小写像の原理 定理 5.1 p.94
𝑹𝑹𝒏𝒏 上で定義された関数𝒈𝒈(𝒙𝒙) が(1) 𝒙𝒙∈𝑹𝑹𝒏𝒏→𝒈𝒈(𝒙𝒙)∈𝑹𝑹𝒏𝒏(2) 𝒙𝒙, 𝒚𝒚∈𝑹𝑹^𝒏𝒏→||𝒈𝒈(𝒙𝒙)−𝒈𝒈(𝒚𝒚)||≤𝑳𝑳||𝒙𝒙−𝒚𝒚||∈𝑹𝑹^𝒏𝒏(3) 𝟎𝟎≤𝑳𝑳<𝟏𝟏を満たすとき，方程式
𝒙𝒙=𝒈𝒈(𝒙𝒙)の解𝒙𝒙=𝜶𝜶は𝑹𝑹𝒏𝒏においてただ 1 つ存在する．そして，𝒙𝒙(𝟎𝟎) を「𝒙𝒙(𝒌𝒌+𝟏𝟏)=𝒈𝒈(𝒙𝒙(𝒌𝒌))によって，生成される列〖{𝒙𝒙〗^((𝒌𝒌))}は𝒌𝒌→∞のとき𝜶𝜶に「1」とする反復に」する．

スライド５ 縮小写像の原理 定理 5.2 p.94
あるノルム||・||について|(|𝐌𝐌|)|<𝟏𝟏ならば反復法
𝒙𝒙^((𝒌𝒌+𝟏𝟏))=𝝓𝝓(𝒙𝒙^((𝒌𝒌) ) )=𝑴𝑴𝑴𝑴^((𝒌𝒌) )+𝑵𝑵𝑵𝑵によって作られる列〖{𝒙𝒙〗^((𝒌𝒌))}は(5.1)の解𝒙𝒙に収束する(証明)

任意の𝒙𝒙, 𝒚𝒚∈𝑹𝑹^𝒏𝒏に対して，|(|𝝓𝝓(𝒙𝒙)−𝝓𝝓(𝒚𝒚)|)|=|(|𝑴𝑴𝑴𝑴−𝑴𝑴𝑴𝑴|)|≤|(|𝑴𝑴|)|∙||𝒙𝒙−𝒚𝒚||(≤||𝒙𝒙−𝒚𝒚|)|が成り立つので，𝝓𝝓(𝒙𝒙)は定理 5.1 の仮定を満たす．よって， 〖{𝒙𝒙〗^((𝒌𝒌))}は解𝒙𝒙に収束する．
スライド 6 行列 A の分解

スライド 7 ヤコビ法 p.95

スライド 8 ヤコビ法：プログラム p.96
#include <stdio.h>#include <stdlib.h>#include <math.h>#define N
10/* N 元方程式 */#define EPS pow(10.0, -8.0)/* epsilon の設定 */#define KMAX 100/* 最大反復回数 */int double_comp(const void *s1 , const void *s2 );/* 最大値ノルムの計算 a[0,...,n-1] */double vector_norm_max(double a[N]);void input_matrix( double a[N][N], char c, FILE *fin, FILE *fout);
void input_vector( double b[N], char c, FILE *fin, FILE *fout);2

/* ヤコビ法 */void jacobi_lin(double a[N][N], double b[N], double x[N]);input_sp.dat
int main(void)
5.0 2.0 0.0 0.0 0.0{
FILE *fin, *fout;
2.0 5.0 2.0 0.0 0.0
double a[N][N], b[N], x[N];0.0 2.0 5.0 2.0 0.0
int i;/* ファイルのオープン */0.0 0.0 2.0 5.0 2.0
if ( (fin = fopen( "input_sp.dat", "r")) == NULL )0.0 0.0 0.0 2.0 5.0{
printf("ファイルが見つかりません : input_sp.dat ¥n");0.0 0.0 0.0 0.0 2.0
exit(1);
0.0 0.0 0.0 0.0 0.0}
if( (fout = fopen( "output_sp.dat", "w")) == NULL )0.0 0.0 0.0 0.0 0.0{
0.0 0.0 0.0 0.0 0.0
printf("ファイルが作成できません : output_sp.dat ¥n");exit(1);
0.0 0.0 0.0 0.0 0.0}
3.0 1.0 4.0 0.0 5.0
input_matrix( a, 'A', fin, fout ); /* 行列 A の入出力 */input_vector( b, 'b', fin, fout ); /* ベクトル b の入出力 */1.0 1.0 1.0 1.0 1.0
input_vector( x, 'x', fin, fout ); /* 初期ベクトル x0 の入出力*/
jacobi_lin( a, b, x );/* ヤコビ法 *//* 結果の出力 */fprintf( fout, "Ax=b の解は次の通りです¥n");for( i = 0 ; i < N ; i++){fprintf(fout, "%f¥n", x[i]);}fclose(fin); fclose(fout); /* ファイルのクローズ */return 0;}
void jacobi_lin(double a[N][N], double b[N], double x[N])output_sp.dat{
行列 A は次の通りです
double eps, xn[N];
5.00
2.00
0.00
int i, j, k=0;
2.00
5.00
2.00
/* x <- x_k, xn <- x_{k+1} */（省略）
do{for(i = 0; i < N; i++)ベクトル b は次の通りです{
3.00
xn[i] = b[i];1.00
for ( j = 0; j < N; j++)
4.00
{
0.00
xn[i] -= a[i][j]*x[j];
5.00
}
-1.00
xn[i] += a[i][i]*x[i];/* 余分に引いた分を加える */6.00
xn[i] /= a[i][i];
-2.00
}
7.00
for( i = 0; i < N; i++) x[i] = xn[i]-x[i];-15.00eps = vector_norm_max(x);ベクトル x は次の通りですfor( i = 0; i < N; i++) x[i] = xn[i]; /* 値を更新 */
1.00
k++;
1.00
}while(eps > EPS && k < KMAX);
1.00
if ( k == KMAX ){1.00
printf("答えが見つかりませんでした¥n");
1.00
exit(1);
1.00
}
1.00
else{
1.00
printf("反復回数は%d 回です¥n", k);
1.00
3

1.00
Ax=b の解は次の通りです
1.000000-1.000000
2.000000

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

2.0

0.0

0.0

0.0

0.0

5.0

2.00.0

0.0

0.0

2.0

5.0

2.0

0.0

0.0

0.0

2.0

5.0

2.0

0.0

0.0

0.0

2.0

5.0

2.0

0.0

0.0

0.0

2.0

5.0

-1.0

6.0

-2.0

7.0

-15.0

1.0

1.0

1.0

1.0

1.0

0.00
0.00

…
…

0.00
0.00

}

/* 反復回数を画面に表示 */}/* 最大値ノルムの計算 a[0,...,n-1] */double vector_norm_max( double a[N]){int i;for ( i = 0 ; i < N; i++) a[i] = fabs(a[i]);/* 並べ換え：ソートする関数 qsort() */qsort(a, N, sizeof(a[0]), double_comp);/* for(i=m;i<=n;i++) printf("a[%d]=%f¥n",i,a[i]);*/
return a[0];}int double_comp( const void *s1 , const void *s2 ){const double a1 = *((double *)s1); /* (double *)へキャスト */const double a2 = *((double *)s2); /* (double *)へキャスト */if( a1 < a2 ){return -1;}
else if( a1 == a2 ){return 0;}
else{return 1;}}/* 反復回数 66 回 */スライド 9 ガウス・ザイデル法

p.100

反復法
反復回数 34 回
スライド 10 ガウスザイデル法：プログラム p.101/* ガウス・ザイデル法 */double* gauss_seidel(double** a, double* b, double* x){
double eps, * xo, s, t;
int i, j, k = 0;xo = dvector(1, N); /* xo[1...N] */do{/* xo <- x_k, x <- x_{k+1} */for (i = 1; i <= N; i++) xo[i] = x[i]; /* x_k に x_(k+1)を代入 *//* i=1 の処理 */t = 0.0;4}for (j = 2; j <= N; j++) t += a[1][j] * xo[j];x[1] = (b[1] - t) / a[1][1];/* i=2,3, . . .N の処理 */for (i = 2; i <= N; i++){s = 0.0; t = 0.0;
for (j = 1; j < i; j++) s += a[i][j] * x[j];/* i-1 列までの和 */for (j = i+1; j <= N; j++) t += a[i][j] * x[j];/* i+1 列以降の和 */x[i] = (b[i] - s - t) / a[i][i];}for (i = 1; i <= N; i++) xo[i] = xo[i] - x[i];eps = vector_norm_max(xo, 1, N);k++;} while (eps > EPS && k < KMAX);
free_dvector(xo, 1); /* */if (k == KMAX){
printf("答えはがみつかりませんでした¥n");exit(1);}
else{
printf("反復回数は%d です¥n", k); /* 反復回数を画面に表示 */return x;}

スライド 11 データ入出力 p.19
プログラム 2.3
#include <stdio.h>#include <stdlib.h>#define N 4
void input_matrix(double **a,char c,FILE* fin,FILE* fout);
void input_vector(double *b,char c,FILE* fin,FILE* fout);
double **dmatrix(double **a,int nr1,int nr2,int nl1,int nl2);
Void free_dmatrix(double **a,int nr1,int nr2,int nl1,int nl2);
double *dvector(double *a, int i);
Void free_dvector(double *a, int i);
int main(void){FILE *fin,*fout;
double **a,*b;
a = dmatrix(1, N, 1, N); /* A[1…N][1…N] */b = dvector(1,N); /* b[1…N] */if((fin=fopen("input.dat","r"))==NULL){
printf("ファイルは見つかりません:input.dat ¥n");exit(1);}if((fout=fopen("output.dat","w"))==NULL){
printf("ファイルが作成できません:output.dat ¥n");exit(1);}
input_matrix(a,’A’,fin,fout); /* 行列 A の入出力 */input_vector(b,’b’,fin,fout); /* ベクトル b の入出力 */fclose(fin); fclose(fout);}
void input_matrix(double **a,char c,FILE* fin,FILE* fout){
int i,j;fprintf(fout,"行列%c は次の通りです¥n",c);5

for(i=0;i<N;++i){
for(j=0;j<N;++j){
fscanf(fin,"%lf",&(a[i][j]));fprintf(fout,"%5.2f¥t",a[i][j]);}fprintf(fout,"¥n");}}void input_vector(double *b,char c,FILE* fin,FILE* fout){
int i;
fprintf(fout,"ベクトル%c は次の通りです¥n",c);for(i=0;i<N;++i){
fscanf(fin,"%lf",&(b[i]));fprintf(fout,"%5.2f¥t",b[i]);fprintf(fout,"¥n");}}6