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

連立 1 次方程式の直接解法•「•

部分「」消去法」選択付きガウス消去法

スライド２ 連立 1 次方程式 p.38
•連立 1 次方程式の重要性•「•

非線形問題を「」化して解く．

•

複雑な構造物，「」な要素に分解•

各要素に対して「•「•

ただし，変数の数は「」•

コンピュータにより「」に解く．

」の問題は基本的には解けない．

」方程式を立てる．
」の原理により統合

スライド３ ガウス消去法 p.38
連立 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
ガウス消去法は，「」をコンピュータに適した方法で行う．

スライド４ ガウス消去法 p.38

(1) 第 1 列の 2 行目と３行目を 0 にする．
「第 1 行×(－5/3)+第 2 行目」「第 1 行×(－4/3)+第 3 行目」(2) 第 2 列の 3 行目を 0 にする．
「第 2 行×1+第 3 行」
(1), (2) : 「」消去

(3) x3, x2, x1 の順に代入して答えを求める．
(3) : 「」消去
1

スライド５ ガウスの消去法：ｎ元（「p.39」）スライド６ ガウス消去法：ｎ元p.40

スライド７ ガウス消去法：アルゴリズム p.41
行列 A とベクトル b の入力/* 前進消去 *//* 後退代入 */bを/* 答えは b に上書き*/

出力2スライド８ ガウス消去法：プログラム p.41
#include <stdio.h>#include <stdlib.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 simple_gauss(double a[N][N],double b[N]);int main(void){
input.dat の内容
FILE *fin, *fout;
1.0
2.0
1.0
1.0
double a[N][N], b[N];
4.0
5.0
-2.0
4.0
int i;
if((fin=fopen("input.dat","r"))==NULL) exit(1);
4.0
3.0
-3.0
1.0if((fout=fopen("output.dat","w"))==NULL) exit(1);
2.0
1.0
1.0
3.0
input_matrix(a,’A’,fin,fout); input_vector(b,’b’,fin,fout);
-1.0
-7.0
-12.0
2.0
simple_gauss(a,b);
fprintf(fout,"Ax=b の解は次の通りです¥n");for(i=0;i<N;i++){ fprintf(fout,"%f¥n",b[i]);}output.dat の内容
fclose(fin); fclose(fout);
return 0;
行列 A は次の通りです}
1.00
2.00
1.00
1.00
void simple_gauss(double a[N][N],double b[N]){
4.00
5.00
-2.00
4.00
int i,j,k;double alpha, tmp;
4.00
3.00
-3.00
1.00
for(k=0;k<N-1;k++){ /* 前進消去 */
2.00
1.00
1.00
3.00
for(i=k+1;i<N;i++){alpha=-a[i][k]/a[k][k];ベクトル b は次の通りですfor(j=k+1;j<N;j++){-1.00
a[i][j]=a[i][j]+alpha*a[k][j];-7.00}
b[i]=b[i]+alpha*b[k];-12.00}2.00}
Ax=b の解は次の通りですb[N-1]=b[N-1]/a[N-1][N-1]; /* 後退代入 */for(k=N-2;k>=0;k--){-2.000000
tmp=b[k];-1.000000
for(j=k+1;j<N;j++){tmp=tmp-a[k][j]*b[j];1.000000}
2.000000b[k]=tmp/a[k][k];}}3スライド９ 部分ピボット付きガウス消去法 p.44

(1) 第 1 列の絶対最大要素:第 3 行．第 1 行と第 3 行を交換
(2) 第 1 列の第 2 行，第 3 行を 0 にする．
(3) 第 2 列以下で第 2 列の絶対最大要素:第 3 行．第 2 行と第 3 行を交換
(4) 第 2 列の第 3 行を 0 にする．
スライド１０ 部分ピボット選択のアルゴリズム p.45
n 元連立 1 次方程式､すでに消去作業を k-1 回行っていると仮定
1) k 回目の消去作業において|aik|(i=k,k+1,…,n-1)のうちで最大のものを調べて，その行番号を ipとする．
2) ip≠k ならば，k 行目と ip 行目を入れ換えて消去作業を行う．
アルゴリズム続き

スライド１１ データ入出力 p.19 (「C 言語入門」p.282)
プログラム 2.3 改#include <stdio.h>#include <stdlib.h>#define N 4
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);

4int main(void){
FILE *fin,*fout;
double a[N][N],b[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.dat の内容
1.0 2.0 1.0 1.0
4.0 5.0 -2.0 4.0
4.0 3.0 -3.0 1.0}
input_matrix(a,’A’,fin,fout); /* 行列 A の入出力 */input_vector(b,’b’,fin,fout); /* ベクトル b の入出力 */2.0 1.0 1.0 3.0-1.0 -7.0 -12.0 2.0

fclose(fin); fclose(fout);}
void input_matrix(double a[N][N],char c,FILE* fin,FILE* fout){
int i,j;fprintf(fout,”行列%c は次の通りです¥n”,c);for(i=0;i<N;++i){

output.dat の内容
行列 A は次の通りです．

for(j=0;j<N;++j){1.0 2.0 1.0 1.0

fscanf(fin,”%lf”,&(a[i][j]));4.0 5.0 -2.0 4.0

fprintf(fout,”%5.2f¥t”,a[i][j]);4.0 3.0 -3.0 1.0}

2.0 1.0 1.0 3.0

fprintf(fout,”¥n”);

ベクトル b は次の通りです．

}

-1.0}
void input_vector(double b[N],char c,FILE* fin,FILE* fout){
int i,j;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”);
}
}

5-7.0
-12.0
2.0