数値解析
2023年度前期 第9週 [6月15日]

静岡大学大学院
工学研究科機械工学専攻
ロボット・計測情報分野
創造科学技術大学院情報科学専攻
三浦 憲二郎講義アウトライン [6月15日]•連立1次方程式の反復解法•密行列と疎行列•反復法の原理：縮小写像の原理•ヤコビ法•ガウス・ザイデル法•データ入出力

2重のポインターの用法

連立一次方程式の反復解法 p.93
•連立1次方程式
𝑨𝑨𝑨𝑨 = 𝒃𝒃•密(dense)行列と疎（sparse）行列
密行列：行列の要素のすべて，あるいはほとんどゼロでない行列
疎行列：行列の要素のほとんどがゼロである行列•疎行列•反復解法を用いるほうが速い．
•ヤコビ法とガウスザイデル法

反復法の原理 p.93
連立1次方程式
𝑨𝑨𝑨𝑨 = 𝒃𝒃 (5.1)

を同値な方程式

𝒙𝒙 = 𝝓𝝓 𝒙𝒙 = 𝑴𝑴𝑴𝑴 + 𝑵𝑵𝑵𝑵に変形し，反復法により解を求める．𝑴𝑴を反復行列という．

縮小写像の原理

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

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

任意の𝒙𝒙, 𝒚𝒚 ∈ 𝑹𝑹𝒏𝒏 に対して，𝝓𝝓 𝒙𝒙 − 𝝓𝝓 𝒚𝒚= 𝑴𝑴𝑴𝑴 − 𝑴𝑴𝑴𝑴 ≤ 𝑴𝑴 � | 𝒙𝒙 − 𝒚𝒚 ≤ | 𝒙𝒙 − 𝒚𝒚が成り立つので，𝝓𝝓 𝒙𝒙 は定理5.1の仮定を満たす．よって， {𝒙𝒙(𝒌𝒌) }は解𝒙𝒙に収束する．

行列Aの分解 p.95
𝑨𝑨 = 𝑫𝑫 + 𝑳𝑳 + 𝑼𝑼
𝒂𝒂𝟏𝟏𝟏𝟏
𝑫𝑫 = ⋮𝟎𝟎
𝑳𝑳 = ⋮𝒂𝒂𝒏𝒏𝒏𝒏

𝟎𝟎
𝑼𝑼 = ⋮𝟎𝟎⋯⋱⋮⋯ 𝒂𝒂𝒏𝒏𝒏𝒏
⋯ 𝟎𝟎
⋱ ⋮
⋯ 𝟎𝟎

⋯ 𝒂𝒂𝟏𝟏𝟏𝟏⋱⋮⋯ 𝟎𝟎

ヤコビ法 p.95
𝑴𝑴 = −𝑫𝑫−𝟏𝟏 𝑳𝑳 + 𝑼𝑼 , 𝑵𝑵 = 𝑫𝑫−𝟏𝟏とおいた反復法

𝒙𝒙(𝒌𝒌+𝟏𝟏) = 𝑫𝑫−𝟏𝟏 {𝒃𝒃 − (𝑳𝑳 + 𝑼𝑼)𝒙𝒙(𝒌𝒌) }ここで， 𝑫𝑫−𝟏𝟏 の存在を仮定している．
アルゴリズムInput 𝑨𝑨, 𝒃𝒃, 𝒙𝒙(𝟎𝟎) , 𝜺𝜺, 𝒌𝒌𝒌𝒌𝒌𝒌𝒌𝒌
𝒌𝒌 ← 𝟎𝟎

DoFor 𝒊𝒊 = 𝟏𝟏, 𝟐𝟐, … , 𝒏𝒏(𝒌𝒌+𝟏𝟏)𝒙𝒙𝒊𝒊

end for←𝟏𝟏(𝒌𝒌)(𝒃𝒃𝒊𝒊 − ∑𝒏𝒏𝒋𝒋≠𝒊𝒊 𝒂𝒂𝒊𝒊𝒊𝒊 𝒙𝒙𝒋𝒋 )𝒂𝒂𝟏𝟏𝟏𝟏ヤコビ法 p.95
アルゴリズムInput 𝑨𝑨, 𝒃𝒃, 𝒙𝒙(𝟎𝟎) , 𝜺𝜺, 𝒌𝒌𝒌𝒌𝒌𝒌𝒌𝒌
𝒌𝒌 ← 𝟎𝟎

DoFor 𝒊𝒊 = 𝟏𝟏, 𝟐𝟐, … , 𝒏𝒏𝒙𝒙(𝒌𝒌+𝟏𝟏) ←end for𝒌𝒌 ← 𝒌𝒌 + 𝟏𝟏

𝟏𝟏(𝒌𝒌)(𝒃𝒃𝒊𝒊 − ∑𝒏𝒏𝒋𝒋≠𝒊𝒊 𝒂𝒂𝒊𝒊𝒊𝒊 𝒙𝒙𝒋𝒋 )𝒂𝒂𝟏𝟏𝟏𝟏

while ( 𝒙𝒙 𝒌𝒌+𝟏𝟏 − 𝒙𝒙 𝒌𝒌+𝟏𝟏> 𝜺𝜺 and 𝒌𝒌 < 𝒌𝒌𝒎𝒎𝒎𝒎𝒎𝒎)

If 𝒌𝒌 = 𝒌𝒌𝒎𝒎𝒎𝒎𝒎𝒎 “失敗” else Output 𝒙𝒙 𝒌𝒌+𝟏𝟏 end if

ヤコビ法：プログラム その１#include <stdio.h>#include <stdlib.h>#include <math.h>#define N#define EPS#define KMAX

10/* N元方程式 */pow(10.0, -8.0)/* epsilonの設定 */100/* 行列の入力 */void input_matrix(double **a, char c, FILE* fin, FILE* fout);/* ベクトルの入力 */void input_vector(double *b, char c, FILE* fin, FILE* fout);/* 行列の領域確保 */double **dmatrix(int nr1, int nr2, int nl1, int nl2);/* 行列の領域解法 */void free_dmatrix(double** a, int nr1, int nr2, int nl1, intnl2);/* ベクトル領域の確保 */double* dvector(int i, int j);/* ベクトルの領域解法 */void free_dvector(double* a, int i);/* 比較関数 */int double_comp(const void *s1, const void *s2);/* 最大値ノルムの計算 a[m...n] */double vector_norm_max(double *a, int m, int n);/* ヤコビ法 */double* jacobi_lin(double** a, double* b, double* x);ヤコビ法：プログラム その２
int main(void){
FILE* fin, * fout;
double** a, * b, * x;
int i;/* 行列およびベクトルの領域確保 */a = dmatrix(1, N, 1, N); /* 行列 a[l. . .N] [1. . .N] */b = dvector(1, N); /* b[l...N] */x = dvector(1, N); /* x[l...N] *//* ファイルのオープン */if ((fin = fopen("input_sp.dat", "r")) == NULL){
printf("ファイルが見つかりません : input_sp.dat ¥n");exit(1);}if ((fout = fopen("output_sp.dat", "w")) == NULL){
printf("ファイルが作成できません : output_sp.dat ¥n");exit(1);}
input_matrix(a, 'A', fin, fout); /* 行列Aの入出力 */input_vector(b, 'b', fin, fout); /* ベクトルbの入出力 */input_vector(x, 'x', fin, fout); /* 初期ベクトルx0の入出力 */x = jacobi_lin(a, b, x);/* 結果の出力 */fprintf(fout, "Ax=b の解は次の通りです¥n");for (i = 1; i <= N; i++){fprintf(fout, "%f¥n", x[i]);}fclose(fin); fclose(fout); /* ファイルのクローズ *//* 領域の解法 */free_dmatrix(a, 1, N, 1, N); free_dvector(b, 1);
free_dvector(x, 1);
return 0;}ヤコビ法：プログラム その３/* ヤコビ法 */double* jacobi_lin(double** a, double* b, double* x){
double eps, *xn;
int i, j, k = 0;xn = dvector(1, N); /* xn[1...N] *//* x <- x_k, xn <- x_{k+l} */do{for (i = 1; i <= N; i++){xn[i] = b[i];for (j = 1; j <= N; j++){xn[i] -= a[i][j] * x[j];}xn[i] += a[i][i] * x[i];xn[i] /= a[i][i];}for (i = 1; i <= N; i++) x[i] = xn[i] - x[i];eps = vector_norm_max(x, 1, N); /* 最大値ノルムの計算 */for (i = 1; i <= N; i++) x[i] = xn[i]; /* 値を更新 */k++;} while (eps > EPS && k < KMAX);
free_dvector(xn, 1); /* 領域の解法 */if (k == KMAX){
printf("答えが見つかりませんでした¥n");exit(1);}
else{
printf("反復回数は%dです¥n", k); /* 反復回数を画面に表示 */return x;}}ガウス・ザイデル法 p.100
とおいた反復法

𝑴𝑴 = − 𝑫𝑫 + 𝑳𝑳 −𝟏𝟏 𝑼𝑼, 𝑵𝑵 = 𝑫𝑫 + 𝑳𝑳 −𝟏𝟏𝒙𝒙(𝒌𝒌+𝟏𝟏) = 𝑫𝑫−𝟏𝟏 {𝒃𝒃 − 𝑳𝑳𝒙𝒙 𝒌𝒌+𝟏𝟏 − 𝑼𝑼𝒙𝒙 𝒌𝒌 }ここで， 𝑫𝑫−𝟏𝟏 の存在を仮定している．
アルゴリズムInput 𝑨𝑨, 𝒃𝒃, 𝒙𝒙(𝟎𝟎) , 𝜺𝜺, 𝒌𝒌𝒌𝒌𝒌𝒌𝒌𝒌
𝒌𝒌 ← 𝟎𝟎

DoFor 𝒊𝒊 = 𝟏𝟏, 𝟐𝟐, … , 𝒏𝒏(𝒌𝒌+𝟏𝟏)𝒙𝒙𝒊𝒊

end for←𝟏𝟏𝒌𝒌+𝟏𝟏(𝒃𝒃𝒊𝒊 − ∑𝒊𝒊−𝟏𝟏𝒂𝒂𝒙𝒙
𝒋𝒋=𝟏𝟏 𝒊𝒊𝒊𝒊 𝒋𝒋
𝒂𝒂𝟏𝟏𝟏𝟏(𝒌𝒌)− ∑𝒏𝒏𝒋𝒋≠𝒊𝒊 𝒂𝒂𝒊𝒊𝒊𝒊 𝒙𝒙𝒋𝒋 )ガウスザイデル法：プログラム 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;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;}}

データ入出力 p.19 (「明快入門C」p.161)
プログラム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] */データ入出力 p.19 (「明快入門C」p.161)
プログラム2.3 続き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);

データ入出力 p.19 (「明快入門C」p.161)
プログラム2.3 続き
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);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");データ入出力 p.19 (「明快入門C」p.161)
プログラム2.3 続き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");}input.datの内容
1.0

2.0

4.03.0 -3.0

4.0

2.0

1.0

1.0

5.0 -2.0

4.0

1.0

3.0

1.0

1.0

-1.0 -7.0 -12.0 2.0

データ入出力 p.19 (「明快入門C」p.161)
プログラム2.3改 続き
output.datの内容
行列Aは次の通りです．
1.0

4.0

4.0

2.0

2.0

1.0

5.0 -2.0

3.0 -3.0

1.0

1.0

1.0

4.0

1.0

3.0

ベクトルbは次の通りです．
-1.0-7.0-12.0
2.0

まとめ
• 連立1次方程式の反復解法•密行列と疎行列•反復法の原理：縮小写像の原理•ヤコビ法•ガウス・ザイデル法•データ入出力

2重のポインターの用法