1. No category

数値シミュレーション 第 4 回 教材
担 当 神保 雅人
2014 年 10 月 17 日（金）実施
連立 1 次方程式の数値解
ガウスの消去法と後退代入
ここでは，簡単のため，変数 x1，x2，x3 についての 3 元連立 1 次方程式を考える。
a x +a x +a
a x +a x +a
a x +a x +a
x =b
x =b
x =b
(1)
(2)
(3)
(1)式の両辺を a11 で割ると，
x +
a
a
x +
a
a
x =
b
a
(4)
(4)式に a21 を掛けたものを(2)式から引いて a21 を消去する[(5)式]。また，(4)式に a31 を掛けた
ものを(3)式から引いて a31 を消去する[(6)式]。
x
a
a
+
a a
a
a a
−
a
−
a
a
x
x + a
a
a
a a
−
a
a a
−
a
+
x + a
x =
x +a
x +a
x +a
x =b
x =b
x =b
(4)
b a
a
b a
x =b −
a
x =b −
この結果を次の様に書き直す。
x +a
a
a
b
a
(5)
(6)
(1' )
(2' )
(3' )
(2' )式の両辺を a22' で割り[(7)式]，それに a32' を掛けたものを(3' )式から引いて a32' を消去す
る[(8)式]。
x +a
a
x +a
a
x +
a
a a
−
a
x =b
b
x =
a
x =b −
(1' )
(7)
b a
a
(8)
更に，この結果を次の様に書き直す。
x +a
x +a
x +a
a
x =b
x =b
x =b
(1' )
(2'' )
(3'' )
－ 1 －
数値シミュレーション 第 4 回 教材
担 当 神保 雅人
この様に aij (i > j)を消去する手続きをガウスの消去法(Gaussian elimination)と呼ぶ。
ここで，(3'' )式から x3 について解が求まる[(9)式]。
x =
b
a
(9)
これを(2'' )式に代入すれば，x2 について解が求まる[(10)式]。
x =b
−a
x
(10)
更に，x2 及び x3 の解を(1' )式に代入すれば，x1 について解が求まる[(11)式]。
x = b − (a
x +a
x )
(11)
この様に，得られた解を一行前の方程式に代入して，順に解を求めていく手続きを後退代入
(back substitution) と呼ぶ。
連立 1 次方程式は，ガウスの消去法を適用した後，後退代入を適用することによって，数値的
に解くことが出来る。但し，ガウスの消去法は，計算途中で対角成分が 0 に近くなると適用出来
なくなるので，列の入れ替えが必要となる。
一般の n 元連立 1 次方程式の扱いや，他の数値解法については【文献 4-1】を，様々な数値解
法の計算回数の比較や C 言語プログラムについては【文献 4-2】を参照のこと。
【文献 4-1】森正武，
『数値解析 第 2 版』
，共立出版，2002 年 [第 1 版は 1973 年出版]
【文献 4-2】W. H. Press et al，丹慶勝市 他 訳，
『Numerical Recipes in C [日本語版]』
，
技術評論社，1993 年
例題 1
次のプログラムはガウスの消去法と後退代入とを用いて，3 元連立 1 次方程式を数値的に解く
ものである。これを入力，ビルドして，実行せよ。ここで，クラス名は Sample4_1，ソースファ
イル名は Sample4_1.java とする。
なお，ここで解く方程式は次のものであり，プログラムの中の配列の添え字は 0 から始めるた
めに，一つずつずらしてある。
【 c[i-1][j-1] = aij (i,j=1,2,3), c[i-1][3] = bi (i=1,2,3) 】
3x + x − x = 5
2x − 2x + x = 6
4x + 3x − 2x = 7
public class Sample4_1 {
public static void main(String[] args) {
final int n = 3;
double[][] c = {{3.0, 1.0, -1.0, 5.0}, {2.0, -2.0, 1.0, 6.0}, {4.0, 3.0, -2.0,
7.0}};
double[] x = new double[n];
Gaussel(c, n);
－ 2 －
数値シミュレーション 第 4 回 教材
担 当 神保 雅人
backsub(c, x);
for (int i=0; i<n; i++) {
for (int j=0; j<n+1; j++)
System.out.print(" c_"+(i+1)+(j+1)+"="+c[i][j]);
System.out.println();
}
for (int i=0; i<n; i++)
System.out.println("x_"+(i+1)+"="+x[i]);
}
public static void Gaussel(double[][] c, int n) {
double cii, cki;
for (int i=0; i<n; i++) {
cii = c[i][i];
for (int j=i; j<n+1; j++)
c[i][j] = c[i][j]/cii;
for (int k=i+1; k<n; k++) {
cki = c[k][i];
for (int m=i; m<n+1; m++)
c[k][m] = c[k][m] - c[i][m]*cki;
}
}
}
public static void backsub(double[][] c, double x[]) {
int n = x.length;
for (int i=n-1; i>=0; i--) {
double sum = 0.0;
for (int j=i+1; j<n; j++)
sum = sum + c[i][j]*x[j];
x[i] = c[i][n] - sum;
}
}
}
演習 1
例題 1 の 3 元連立 1 次方程式を数値的に解いた結果が元の方程式を満たしていることを確かめ
るプログラムを作成し，ビルドして，実行せよ。ここで，クラス名は Ex4_1，ソースファイル名
は Ex4_1.java とする。(ヒント：例題のプログラムは x1，x2，x3 の解を表示させているが，これ
らを連立方程式の第 1 式，2 式，3 式の左辺に代入した結果を表示して，右辺の値となるかどうか
を確かめればよい)
－ 3 －
数値シミュレーション 第 4 回 教材
担 当 神保 雅人
今回の提出物
1) Ex4_1.java
2) 演習 1 の出力結果(Eclipse のコンソールから文字列をコピー)
－ 4 －