この記事は弊学のプログラミングの授業の内容を自分なりにまとめたものです。
間違い・不適切な表現があればこのTwitterアカウントにDM送ってください。
上三角形型の方程式の解法
人間が解く
まずは次のような連立方程式を解く
\begin{equation}
\left.
\begin{aligned}
2x-y+0z-3u=1 \\
2y+3z+7u=0 \\
z+9u=6 \\
5u=-3 \\
\end{aligned}
\right\}
\
\end{equation}
このような連立方程式は比較的容易に解くことが出来る事は一般に知られています。
まず行列の対角成分を1にする為にそれぞれの式の対角成分となる変数の係数で割っていきます。
これを行うと以下のようになります。
\begin{equation}
\left.
\begin{aligned}
x-0.5y+0z-1.5u=0.5 \\
y+1.5z+3.5u=0 \\
z+9u=6 \\
u=-0.6 \\
\end{aligned}
\right\}
\
\end{equation}
ここまでたどり着いたなら下から各変数の値を求め、その変数をそれより上の式に代入する値を代入するという行為を繰り返すだけで解を求めることができます。
プログラムが解く
これをプログラムで実装したいのですが拡大係数行列で入力できるようにしたいと思います。
上の連立方程式を拡大係数行列で表現すると以下のようになります。
\begin{pmatrix}
1 & 0.5 & 0 & 1.5 & 0.5 \\
& 1 & 1.5 & 3.5 & 0\\
& & 1 & 9 & 6\\
& & & 1 & -0.6
\end{pmatrix}
拡大係数行列に関しては今回の話の本質とは脱線してしまうので省かせていただきます。
で、以下のコードが上三角形行列を解くソースコードです。
(あとで以上ずつ解説します)
#include<iostream>
#include<algorithm>
#define ZERO 1e-8
using namespace std;
/*問題要件:
nxnの上三角形行列が与えられる。
それを解くプログラムを作成せよ。
入力
n
a_11 a_12 a_13 ... a_1n-1 a_1n b1
a_22 a_23 ... a_2n-1 a_2n b2
a_33 ... a_3n-1 a_3n b3
...
a_n-1n-1 a_n-1n bn-1
a_nn bn
出力
x_1 x_2 ... x_n-1 x_n
*/
int main(){
int n;
cin >> n;
long double a[100][100];
for(int i=0;i<n;i++){
for(int j=i;j<=n;j++){
cin >> a[i][j];
}
}
for(int i=0;i<n;i++){
int beginNum = a[i][i];
if(!beginNum){
cout << "can not calc" << endl;
return 0;
}
for(int j=i;j<=n;j++)a[i][j] /= beginNum;
}
for(int i=n-1;i>=0;i--){
for(int j=i+1;j<n;j++)a[i][n] -= a[j][n]*a[i][j];
}
for(int i=0;i<n;i++)cout << a[i][n] << endl;
}
解説
それぞれの箇所を簡単に解説します。
//入力部
int n;
cin >> n;
long double a[100][100];
for(int i=0;i<n;i++){
for(int j=i;j<=n;j++){
cin >> a[i][j];
}
}
ここは数値を入力する箇所です。
注意点として今回は上三角形行列なのでn行目以降の入力はn列目の要素から開始します。
for(int i=0;i<n;i++){
int beginNum = a[i][i];//各階層の一番最初の数値。
//もし最初の数値が0だった場合は答えが一意に定まらないので処理を中止する。
if(!beginNum){
cout << "can not calc" << endl;
return 0;
}
//i番目の式の右辺全てを最初の数値で割る
for(int j=i;j<=n;j++)a[i][j] /= beginNum;
}
このソースコードのbeginNumは対角成分です。
ここでは行列の対角成分を1にする為にそれぞれの式の対角成分となる変数の係数で割っていきます。
ただし、対角成分が0だった場合は解が一意に定まらないのでプログラムを強制終了しています。
for(int i=n-1;i>=0;i--){
//(i,i)番目の要素を求めるため既に求まっている変数をb[i]から削除します。
//これを下から行うことによって答えが一意に定まります。
//b[i] = i番目の式の値です。
for(int j=i+1;j<n;j++)a[i][n] -= a[j][n]*a[i][j];
}
ここでは変数を下から求めて行っています。
この処理は逆進代入法とか言われたりします。
基本的に人間が行う方法と基本的には変わらないです。
人間がやる場合は左辺を単位行列にするために行列ごと引いたりしますが、今回の問題の要件は変更先はn列目だけで構わないので不要な行列の修正は行なっていません。
for(int i=0;i<n;i++)cout << a[i][n] << endl;
最後の出力部分です。
左辺が単位行列になっているので左辺を上から表示すれば解が求まります。
連立方程式の解法
人間が解く
先ほど上三角形行列の解法は学びました。
では次に一般の連立一次方程式を解く方法を学びましょう。
\begin{equation}
\left.
\begin{aligned}
x+y-z=1 \\
4x-2y+z=-6 \\
9x+3y+z=9 \\
\end{aligned}
\right\}
\
\end{equation}
ガウスの消去法
この方法によって一般の連立方程式を上三角形型の拡大係数行列に変更します。
手順は以下の通りです。
- 未探索の行列の対角成分の絶対値が最大の行の対角成分が1になるように式を割る
- その他の未探索の行列の1の対象となった列の値が0になるように行列を変形させる
- 1に戻る
という非常に単純なものです。
言葉だけだと理解ができなさそうな人はソースコードを見ましょう。
プログラムが解く
#include<iostream>
#include<algorithm>
#include<cmath>
#define ZERO 1e-8
using namespace std;
/*問題要件:
nxnの拡大係数行列が与えられる。
それを解くプログラムを作成せよ。
入力
n
a_11 a_12 a_13 ... a_1n-1 a_1n b1
a_21 a_22 a_23 ... a_2n-1 a_2n b2
a_31 a_32 a_33 ... a_3n-1 a_3n b3
...
a_n-11 a_n-12 a_n-13 ... a_n-1n-1 a_n-1n bn-1
a_n1 a_n2 a_n3 ... a_nn bn
出力
x_1 x_2 ... x_n-1 x_n
*/
int main(){
int n;
cin >> n;
long double a[100][100];
for(int i=0;i<n;i++){
for(int j=0;j<=n;j++){
cin >> a[i][j];
}
}
for(int i=0;i<n;i++){
int maxLine=i;
for(int j=i;j<n;j++) if(abs(a[maxLine][i])<abs(a[j][i])) maxLine=j;
for(int j=i;j<=n;j++) swap(a[i][j],a[maxLine][j]);
long double beginNum = a[i][i];
if(!beginNum){
cout << "can not calc" << endl;
return 0;
}
for(int j=i;j<=n;j++)a[i][j] /= beginNum;
for(int j=i+1;j<n;j++){
long double beginDelete = a[j][i];
for(int k=i;k<=n;k++)a[j][k] -= beginDelete*a[i][k];
}
}
for(int i=n-1;i>=0;i--){
for(int j=i+1;j<n;j++)a[i][n] -= a[j][n]*a[i][j];
}
for(int i=0;i<n;i++)cout << a[i][n] << endl;
}
解説
//入力部
int n;
cin >> n;
long double a[100][100];
for(int i=0;i<n;i++){
for(int j=0;j<=n;j++){
cin >> a[i][j];
}
}
まずは入力部です。
先ほどはj=iから始まっていましたが今回は逆三角形行列ではないので0から始めます。
for(int i=0;i<n;i++){
//第n列目の要素の絶対値が最大の行を第n行と交換
int maxLine=i;
for(int j=i;j<n;j++) if(abs(a[maxLine][i])<abs(a[j][i])) maxLine=j;//探索
for(int j=i;j<=n;j++) swap(a[i][j],a[maxLine][j]);//交換
long double beginNum = a[i][i];//各階層の一番最初の数値。
//もし最初の数値が0だった場合は答えが一意に定まらないので処理を中止する。
if(!beginNum){
cout << "can not calc" << endl;
return 0;
}
//i番目の式の右辺全てを最初の数値で割る
for(int j=i;j<=n;j++)a[i][j] /= beginNum;
//n+1行目以降の行のn行目の要素のn列目の要素を0にするように行列を変形する。
for(int j=i+1;j<n;j++){
long double beginDelete = a[j][i];
for(int k=i;k<=n;k++)a[j][k] -= beginDelete*a[i][k];
}
}
ここでガウスの消去法が行われています。
//第n列目の要素の絶対値が最大の行を第n行と交換
int maxLine=i;
for(int j=i;j<n;j++) if(abs(a[maxLine][i])<abs(a[j][i])) maxLine=j;//探索
for(int j=i;j<=n;j++) swap(a[i][j],a[maxLine][j]);//交換
まずは第n列目の要素の絶対値が最大の行を第n行と交換を行います。
これをやる理由はn行目の対角成分が0な場合でも他の配列でそれがカバーできる可能性があるためそれに対応するためです。
long double beginNum = a[i][i];//各階層の一番最初の数値。
//もし最初の数値が0だった場合は答えが一意に定まらないので処理を中止する。
if(!beginNum){
cout << "can not calc" << endl;
return 0;
}
//i番目の式の右辺全てを最初の数値で割る
for(int j=i;j<=n;j++)a[i][j] /= beginNum;
ここは逆三角行列とほとんど同じ処理です。
これで対角成分を１にするように式全体を割ります。
//n+1行目以降の行のn行目の要素のn列目の要素を0にするように行列を変形する。
for(int j=i+1;j<n;j++){
long double beginDelete = a[j][i];
for(int k=i;k<=n;k++)a[j][k] -= beginDelete*a[i][k];
}
でここでn+1行目以降の行のn行目の要素のn列目の要素を0にするように行列を変形します。
変形の方法はn列目が消えるようにn行目をn列目を消したい行のn列目の成分倍した値をn列目を消したい行から引きます。
これはn行目の値はこの直前の処理で1になっていることが保証されているからです。
この処理を繰り返すことで逆三角形行列が完成します。
これ以降の処理は逆三角形行列と全く同じなので解説は省きます。
これで連立一次方程式の解法は終了です。
次回はこれを連立方程式ではなく逆行列に適用した場合の記事を書きます。