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

連立方程式の反復解法•前回の復習•「」法

スライド２ 連立 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

スライド 9 ガウス・ザイデル法

p.100

反復法
反復回数 34 回

2スライド８ 共役勾配法 p.109
正定値行列とその性質」分解と同様，対称正定値(symmetric positive definite)行列の場合に有効

共役勾配法は「n 次元正方行列 A が x≠0 である任意のベクトル x に対して n 次元の内積
p.68, 定理 3.9 A を任意のベクトル x (x≠0)に対して Ax≠0 となるような n 次正方行列とする．このとき，B=AtA は対称正定値行列である．の代わりにを解く．

スライド９ 共役勾配法 p.109
定理 5.7 A は n 次「」行列とする． このとき,連立 1 次方程式 Ａｘ=b の解 ｘ は関数を最小にする．
（証明） 任意のベクトル h (h≠0)に対してA は対称なのでスライド１０ 共役勾配法 p.109

A は正定値なので，（h, Ah） > 0 でスライド１１ 共役勾配法

p.109

共役の定義
ベクトル pi とベクトル pj が（pi,Apj）=0, i≠j を満たすとき ベクトル pi と pjに関して「」である，あるい はＡに関して「」するという．

A を n 次対称正定値行列とする．連立 1 次方程式 Ax=b の解 →3

は行列Ａｆ(x)を最小にする x を求める．

スライド１１ 共役勾配法（conjugate gradient：ＣＧ法） p.109
したがってでf(xk+1)は最小. 残差 rk+1

探索方向 pk+1, ただし p0 = r0 とする．

βk は pk+1 と pk が共役になるように定める．

p.114

スライド１２ 公約勾配法

定理 5.8 共役勾配法によって計算される残差ベクトル rk および 方向ベクトル pk に対して，i≠j のとき次が成り立つ．
（証明）1) k=1β0 の決め方よりA が対称なのでスライド１３ 共役勾配法

p.114

（証明）
2) k>=1 で成り立つと仮定
3) 式(1)a) i<k4b) i=k

スライド１４ 共役勾配法

p.114

（証明）
4) 式(2)a) i<kb) i=k
実は

スライド１５ 共役勾配法

p.115

定理 5.9 A は正則な n 次対称正定値行列で，b は n 次元べクトルとする．このとき，共役勾配法の反復を n 回実行すると Ax=b の解が得られる．
（証明）
n 個のベクトル r0, r1,…, rn-1 は互いに直交しているので，これらすべてに直交するベクトル rn は 0 ベクトルとなる．
したがって，ｒn＝ｂ－Ａｘ＝0,A が正則ならば，スライド１６ 共役勾配法：プログラムその１/* 1 ノルムの計算 */double vector_norm1(double* a, int m, int n){int i;double norm = 0.0;
for (i = m; i <= n; i++){norm += fabs(a[i]);}return norm;}/* A ベクトル a[m...n]と b[m...n]の内積を計算する */double inner_product(int m, int n, double* a, double* b){
5int i;double s = 0.0;
for (i = m; i <= n; i++) s += a[i] * b[i];return s;}/* 行列 a[1...N][1...]とベクトル b[1...N]との積 c <- Ab */void matrix_vector_product(double** a, double* b, double* c){double wk;
int i, j;for (i = 1; i <= N; i++){wk = 0.0;
for (j = 1; j <= N; j++){ wk += a[i][j] * b[j]; }c[i] = wk;}}
スライド１６ 共役勾配法：プログラムその１/* 共役勾配法 */double* cg(double** a, double* b, double* x){
double eps, * r, * p, * tmp, alpha, beta, work;
int i, k = 0;r = dvector(1, N); /* r[l...N] */p = dvector(1, N); /* p[l...N] */tmp = dvector(1, N); /* tmp[l...N] */matrix_vector_product(a, x, tmp); /* tmp <- A b */for( i = 1; i <= N; i++){p[i] = b[i] - tmp[i]; r[i] = p[i];}do{/* alpha の計算 */matrix_vector_product(a, p, tmp); /* tmp <- A p_k */work = inner_product(1, N, p, tmp); /* work <- (p,Ap_k)*/alpha = inner_product(1, N, p, r) / work;/* x_{k+l} と r_{k+l}の計算 */for (i = 1; i <= N; i++) x[i] = x[i] + alpha * p[i];for (i = 1; i <= N; i++) r[i] = r[i] - alpha * tmp[i];/* 収束判定 */eps = vector_norm1(r, 1, N);
k++; /* 反復回数の更新 */if (eps < EPS) goto OUTPUT;/* beta と p_{k+l}の計算 */beta = -inner_product(1, N, r, tmp) / work;for (i = 1; i <= N; i++) p[i] = r[i] + beta * p[i];}while (k < KMAX);OUTPUT:;/* 領域の解放 */free_dvector(r, 1); free_dvector(p, 1); free_dvector(tmp, 1);
if (k == KMAX){
printf("答えが見つかりませんでした¥n");exit(1);}
else{
printf("反復回数は%d です¥n", k); /* 反復回数を画面に表示 */return x;}} /* 反復回数 10 回 */6