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

静岡大学
創造科学技術大学院情報科学専攻
工学部機械工学科
計測情報講座
三浦 憲二郎講義アウトライン [6月8日]•連立1次方程式の直接解法•復習• ガウス消去法•連立1次方程式の直接解法•LU分解復習：連立1次方程式•連立1次方程式の重要性•非線形の問題は基本的には解けない．
•非線形問題を線形化して解く．
•複雑な構造物，単純な要素に分解•各要素に対して線形方程式を立てる．
•重ね合わせの原理により統合•ただし，変数の数は膨大•コンピュータにより数値的に解く．

復習：ガウス消去法
連立1次方程式の解法
x-y=1
(1)
x+2y=4(2)１．代入法
式(1)より y=x-1, 式(2)に代入
x+2(x-1)=4, したがって3x=6, よってx=2, 式(1)よりy=1２．加減法
式(2)より式(1)を引く：3y=3, したがってy=1, 式(1)よりx=2ガウス消去法は，加減法をコンピュータに適した方法で行う．

復習：ガウス消去法：n元(前進消去)

復習：ガウス消去法：n元(後退代入)

LU分解 p.49
1.

1 2をLU分解せよ。
3 4ヒント
1 21=
3 4𝑎𝑎
2.

01𝑏𝑏0𝑐𝑐𝑑𝑑

2 6 4をLU分解せよ。
5 7 9

LU分解
1.

1 2
1 0=
3 43 1

1 20 −22.

1 0 2 64
２ ６ ４=
2.5 1 0 −8 −1
５ ７ ９

LU分解
部分ピボット選択付きガウス消去法(1) k番方程式とip番方程式を入れ換える(2) k番方程式をαik倍してi番方程式(i=k+1,k+2,…,n-1)に加える

LU分解

Ｕ：上三角行列，L：下三角行列，P：置換行列LU分解（Pivotしない場合）
フロベニウス(Frobenius)行列
𝑓𝑓𝑘𝑘 = (0, ⋯ , 0, 𝑓𝑓𝑘𝑘+1,𝑘𝑘 , ⋯ , 𝑓𝑓𝑛𝑛−1,𝑘𝑘 )𝑇𝑇
ここに数式を入力します。
10
01−𝑓𝑓𝑘𝑘+1,𝑘𝑘𝐹𝐹𝑘𝑘 = 𝐼𝐼 − 𝑓𝑓𝑘𝑘 𝑒𝑒𝑘𝑘𝑇𝑇 =⋮0
1−𝑓𝑓𝑛𝑛−1,𝑘𝑘
1
0
01
𝑓𝑓𝑘𝑘+1,𝑘𝑘
1𝐹𝐹𝑘𝑘−1 = 𝐼𝐼 + 𝑓𝑓𝑘𝑘 𝑒𝑒𝑘𝑘𝑇𝑇 =⋮01
𝑓𝑓𝑛𝑛−1,𝑘𝑘

LU分解（Pivotしない場合）
フロベニウス(Frobenius)行列
ここに数式を入力します。
100
10
𝑓𝑓1,0−1
𝑓𝑓2,1=𝐹𝐹0−1 ⋯ 𝐹𝐹𝑛𝑛−2⋮⋮0
𝑓𝑓𝑛𝑛−1,0 𝑓𝑓𝑛𝑛−1,1 1
1
ガウス消去法 前進消去の第k段の操作（行の交換が不要の場合）
𝐴𝐴(𝑛𝑛−1) = 𝐹𝐹𝑛𝑛−2 𝐴𝐴(𝑛𝑛−2) = 𝜋𝜋𝑟𝑟 2

ただし， 𝑓𝑓𝑘𝑘 = (0, ⋯ , 0, 𝑓𝑓𝑘𝑘+1,𝑘𝑘 , ⋯ , 𝑓𝑓𝑛𝑛−1,𝑘𝑘 )𝑇𝑇 , 𝑓𝑓𝑖𝑖,𝑘𝑘 =(𝑘𝑘)𝑎𝑎𝑖𝑖,𝑘𝑘(𝑘𝑘)𝑎𝑎𝑘𝑘,𝑘𝑘−1 (𝑛𝑛−1)
𝐴𝐴(𝑛𝑛−1) = 𝐹𝐹𝑛𝑛−2 ⋯ 𝐹𝐹0 𝐴𝐴 ⟺ 𝐴𝐴 = 𝐹𝐹0−1 ⋯ 𝐹𝐹𝑛𝑛−1𝐴𝐴= 𝜋𝜋𝑟𝑟 2LU分解
以下の行列を計算せよ
1 0 0 −1−𝑎𝑎 1 0−𝑏𝑏 0 1
1 00 10 −𝑐𝑐0 10 −𝑎𝑎
1 −𝑏𝑏

1−𝑎𝑎𝑎𝑎𝑎𝑎 − 𝑏𝑏01−𝑐𝑐0 01 00 10 −101

LU分解

1−𝑎𝑎𝑎𝑎𝑎𝑎 − 𝑏𝑏

1 00 10 −𝑐𝑐

1−𝑎𝑎−𝑏𝑏0 0 −11= 𝑎𝑎
1 00 1𝑏𝑏

0 10 −𝑎𝑎
1 −𝑏𝑏

0 0 −1
1= 𝑎𝑎
1 0𝑐𝑐 1𝑏𝑏

0 01 00 10 01
1 0 = −𝑎𝑎0 1
𝑎𝑎𝑎𝑎 − 𝑏𝑏

0 01 0𝑐𝑐 10 0 1 0 0
1
1 0 0 1 0 = 𝑎𝑎
0 1 0 𝑐𝑐 1
𝑏𝑏

0 01 0𝑐𝑐 1

LU分解

Ｕ：上三角行列，L：下三角行列，P：置換行列

LU分解

LU分解
値の格納

LU分解アルゴリズム
Input A, 𝜀𝜀for 𝑘𝑘 = 0,1, ⋯ , 𝑛𝑛 − 2
Exchange 2 rows if necessary
𝑃𝑃𝑘𝑘 ← 𝑖𝑖𝑖𝑖
if 𝑖𝑖𝑝𝑝 ≠ 𝑘𝑘 thenfor 𝑗𝑗 = 𝑘𝑘, 𝑘𝑘 + 1, ⋯ , 𝑛𝑛 − 1𝑎𝑎𝑘𝑘𝑘𝑘 ↔ 𝑎𝑎𝑖𝑖𝑖𝑖,𝑗𝑗
end forend iffor 𝑖𝑖 = 𝑘𝑘 + 1, 𝑘𝑘 + 2, ⋯ , 𝑛𝑛 − 1𝑎𝑎𝛼𝛼 ← − 𝑖𝑖𝑖𝑖
𝑎𝑎𝑘𝑘𝑘𝑘𝑎𝑎𝑖𝑖𝑖𝑖 ← 𝛼𝛼
for 𝑗𝑗 = 𝑘𝑘 + 1, 𝑘𝑘 + 2, ⋯ , 𝑛𝑛 − 1𝑎𝑎𝑖𝑖𝑖𝑖 ← 𝑎𝑎𝑖𝑖𝑖𝑖 + 𝑎𝑎𝑘𝑘𝑘𝑘
end forend forend for
Output 𝐴𝐴, 𝑃𝑃LU分解
LU分解を利用してAx=bを解くアルゴリズム
Input A, 𝑏𝑏, 𝑃𝑃
for 𝑘𝑘 = 0,1, ⋯ , 𝑛𝑛 − 2
𝑏𝑏𝑘𝑘 ⟷ 𝑏𝑏𝑃𝑃𝑘𝑘 𝐴𝐴
for 𝑖𝑖 = 𝑘𝑘 + 1, 𝑘𝑘 + 2, ⋯ , 𝑛𝑛 − 1𝑏𝑏𝑖𝑖 ← 𝑏𝑏𝑖𝑖 + 𝑎𝑎𝑖𝑖𝑖𝑖 𝑏𝑏𝑘𝑘
end forend forfor 𝑘𝑘 = 𝑛𝑛 − 1, 𝑛𝑛 − 2, ⋯ , 0𝑏𝑏𝑘𝑘 ←end for
Output 𝑏𝑏𝑏𝑏𝑘𝑘 −∑𝑛𝑛−1
𝑗𝑗=𝑘𝑘+1 𝑎𝑎𝑘𝑘𝑘𝑘 𝑏𝑏𝑗𝑗𝑎𝑎𝑘𝑘𝑘𝑘

LU分解：プログラム#include <stdio.h>#include <stdlib.h>#include <math.h>#define N 4 /* 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);

void lu_decomp(double a[N][N],int p[N]);Void lu_solve(double a[N][N],double b[N],int p[N]);int main(void){FILE *fin, *fout;

double a[N][N], b[N];int I, p[N]; /* p[0…N-2]を利用, p[N-1]は未使用 */if((fin=fopen("input_lu.dat","r"))==NULL) exit(1);if((fout=fopen("output_lu.dat","w"))==NULL) exit(1);input_matrix(a,’A’,fin,fout); input_vector(b,’b’,fin,fout);
lu_decomp(a,p); lu_solve(a,b,p);fprintf(fout,"Ax=bの解は次の通りです¥n");for(i=0;i<N;i++){ fprintf(fout,"%f¥n",b[i]);}fclose(fin); fclose(fout);}

return 0;LU分解：lu_decomp()
void lu_decomp(double a[N][N],int p[N]){int i,j,k,ip;
double alpha, tmp;
double amax, eps=pow(2.0,-50.0); /* eps=2^{-50}とする */for(k=0;k<N-1;k++){amax=fabs(a[k][k]); ip=k; /* ピボットの選択 */for(i=k+1;i<N;i++){if(fabs(a[i][k])>amax){
amax=fabs(a[i][k]); ip=i;}}if(amax<eps) { /* 正則性の判定 */printf("入力した行列は正則ではない!!¥n"); exit(1);}p[k]=ip;
if(ip!=k){for(j=k;j<N;j++){tmp=a[k][j]; a[k][j]=a[ip][j]; a[ip][j]=tmp;}}
for(i=k+1;i<N;i++){ /* 前進消去 */alpha=a[i][k]/a[k][k];a[i][k]=alpha;for(j=k+1;j<N;j++){a[i][j]=a[i][j]-alpha*a[k][j];}}}}LU分解：lu_solve()
void lu_solve(double a[N][N],double b[N],int p[N]){int i,j,k;double tmp;
for(k=0;k<N-1;k++){tmp=b[k]; b[k]=b[p[k]]; b[p[k]]=tmp; /*右辺の行変換*/
for(i=k+1;i<N;i++){ /* 前進代入 */b[i]=b[i]-a[i][k]*b[k];}}}b[N-1]=b[N-1]/a[N-1][N-1];
for(k=N-2;k>=0;k--){ /* 後退代入 */tmp=b[k];for(j=k+1;j<N;j++){tmp=tmp-a[k][j]*b[j];}b[k]=tmp/a[k][k];}まとめ•連立1次方程式の直接解法•LU分解