２０２３年度前期

数値解析

第６回授業ノート
２０２３年５月２５日

スライド１
•

講義アウトライン連立 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」）

(1)

スライド６

ガウス消去法：ｎ元（「スライド７(1)𝑎𝑎𝑖𝑖𝑖𝑖 = 𝑎𝑎𝑖𝑖𝑖𝑖 + 𝛼𝛼𝑖𝑖0 𝑎𝑎0𝑗𝑗 , 𝑏𝑏𝑖𝑖
p.40

= 𝑏𝑏𝑖𝑖 + 𝛼𝛼𝑖𝑖0 𝑏𝑏0 , 𝑖𝑖 = 1,2, ⋯ , 𝑛𝑛 − 1」）

ガウス消去法：アルゴリズム

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.0
if((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.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);
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_matrix(a,’A’,fin,fout); /* 行列 A の入出力 */input_vector(b,’b’,fin,fout); /* ベクトル b の入出力 */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){
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[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”);}}

4

input.dat の内容
1.0

2.0

1.0

1.0

4.0

5.0 -2.0

4.0

4.0

3.0 -3.01.0

2.0

1.0

3.0

1.0

-1.0 -7.0 -12.0 2.0

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

2.0

1.0

1.0

4.0

5.0 -2.0

4.0

4.0

3.0 -3.0

1.0

2.0

1.0

3.0

1.0

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