第67回 プログラミングについて『連立１次方程式を解こう その１』
連立１次方程式といえば中学生の頃に学校で習ったような記憶があります。プログラムの世界でも分野にもよりますが、連立１次方程式を解くことは、決して避けることはできないものの１つです。特に、自由曲線を扱うもの、力学的な解析を行なうもの、電気回路の解析を行なうものなどでは当たり前のように使用されています。
私はどちらかと言えば、電気系の分野のプログラムを書くことが多いのですが、電気回路の解析はこれまでの自分の範疇ではなかったのでもっぱら自由曲線を扱うときぐらいにしか使用していないのが正直なところです。
電気回路を勉強した方なら誰でもが知っているキルヒホッフの法則というものがあって、回路に流れる電流を求めるときに非常に役に立ってくれます。
┌───┐ ┌───┐
┌┤ Ｒ１ ├┬┤ Ｒ２ ├┐
│└───┘│└───┘│
│ ← ┌┴┐ ← │
│ │Ｒ│ │＋
│↓Ｉ１↑│ │↓Ｉ２↑○Ｅ
│ │３│ │ー
│ └┬┘ │
│ → │ → │
└─────┴─────┘
上の図で、抵抗が３本（Ｒ１，Ｒ２，Ｒ３）と電源（Ｅ）の回路に流れる電流を求めたいときにこのキルヒホッフの法則を使用します。この法則の内容はともかくとして、上の図で、Ｒ１ーＲ３のループ部分を流れる電流をＩ１と仮定し、Ｒ２ーＲ３ーＥのループを流れる電流をＩ２と仮定します。そうすると、次のような連立１次方程式を作ることができます。
Ｒ１＊Ｉ１＋Ｒ３＊（Ｉ１ーＩ２） ＝０
Ｒ２＊Ｉ２＋Ｒ３＊（Ｉ２ーＩ１）ーＥ＝０
この連立１次方程式を変形して、
Ｒ１＊Ｉ１＋Ｒ３＊Ｉ１ーＲ３＊Ｉ２＝０
Ｒ２＊Ｉ２＋Ｒ３＊Ｉ２ーＲ３＊Ｉ１＝Ｅ
さらに、
（Ｒ１＋Ｒ３）＊Ｉ１ ーＲ３ ＊Ｉ２＝０
ーＲ３ ＊Ｉ１＋（Ｒ２＋Ｒ３）＊Ｉ２＝Ｅ
となります。人間ならばここでさらに式を変形して、Ｉ１とＩ２を求めるのですが、計算機はこのままだと非常に解を求めるのは面倒ですので、この式を行列表現してみます。
┌ ┐┌ ┐ ┌ ┐
│Ｒ１＋Ｒ３、 ーＲ３││Ｉ１│ │０│
│ ││ │＝│ │
│ ーＲ３、Ｒ２＋Ｒ３││Ｉ２│ │Ｅ│
└ ┘└ ┘ └ ┘
するとこのように表現することができます。実際にはＩ１とＩ２以外には数値が入りますので、後は数値の計算をすることになります。例えば、Ｒ１を１、Ｒ２を２、Ｒ３を３、Ｅを１０としてみると、
┌ ┐┌ ┐ ┌ ┐
│１＋３、 ー３││Ｉ１│ │ ０│
│ ││ │＝│ │
│ ー３、２＋３││Ｉ２│ │１０│
└ ┘└ ┘ └ ┘
となり、
┌ ┐┌ ┐ ┌ ┐
│ ４、ー３││Ｉ１│ │ ０│ ー（式１）
│ ││ │＝│ │
│ー３、 ５││Ｉ２│ │１０│ ー（式２）
└ ┘└ ┘ └ ┘
となります。この左側の２Ｘ２の行列が
┌ ┐
│１、０│
│ │
│０、１│
└ ┘
となるようにこの行列を変形してみます。
┌ ┐┌ ┐ ┌ ┐
│ １１／５、 ０││Ｉ１│ │ ６│ ー（式３）
│ ││ │＝│ │
│ー３ 、 ５││Ｉ２│ │１０│ ー（式２）
└ ┘└ ┘ └ ┘
式３は（式１）ー（式２）＊（ー３）／５を計算したものです。次に式２ー式１＊（ー３）／（１１／５）を式４にします。
┌ ┐┌ ┐ ┌ ┐
│１１／５、０││Ｉ１│ │６ │ ー（式３）
│ ││ │＝│ │
│０ 、５││Ｉ２│ │２００／１１│ ー（式４）
└ ┘└ ┘ └ ┘
次に式３／（１１／５）を式５、式４／５を式６にします。
┌ ┐┌ ┐ ┌ ┐
│１、０││Ｉ１│ │ ３０／１１│ ー（式５）
│ ││ │＝│ │
│０、１││Ｉ２│ │ ４０／１１│ ー（式６）
└ ┘└ ┘ └ ┘
結果として、
Ｉ１＝ ３０／１１
Ｉ２＝ ４０／１１
と求めることができます。
プログラムで上の行列を計算するときには、そのままの形では計算し辛いので、
┌ ┐
│Ｉ１│
│ │
│Ｉ２│
└ ┘
を省略し、
┌ ┐ ┌ ┐
│ ４、ー３│ │ ０│
│ │＝│ │
│ー３、 ５│ │１０│
└ ┘ └ ┘
とします。これでも行列が２つですので、
┌ ┐
│ ４、ー３、 ０│
│ │
│ー３、 ５、１０│
└ ┘
と１つの行列にして計算します。
これまでやってきたことをもう少し一般化するとn元の行列は、
┌ ┐
│Ａ00、Ａ01、...Ａ0n、Ｖ0│
│Ａ10、Ａ11、...Ａ1n、Ｖ1│
:
:
│Ａn0、Ａn1、...Ａnn、Ｖn│
└ ┘
という n*(n+1) の行列にして、Ａ00 から Ａnn までの対角の成分が１それ以外が０になるように計算していきます。要するに単位行列にする訳です。
これまで説明してきた行列の解法はガウス・ジョルダン法という方法で、プログラムも単純に記述することができます。
ここからは、次の４ｘ４の行列を解くことを考えてみましょう。
┌ ┐┌ ┐ ┌ ┐
│2.0, 5.0, 8.0, 3.0││Xa│ │10.0│
│4.0, 2.0, 3.0, 7.0││Xb│＝│25.0│
│8.0, 6.0, 9.0, 4.0││Xc│ │30.0│
│9.0, 4.0, 3.0, 8.0││Xd│ │45.0│
└ ┘└ ┘ └ ┘
先ほど説明したように、この行列を次のようにプログラム内で表現します。
#define DIM 4
double a[DIM][DIM+1]={
2.0, 5.0, 8.0, 3.0, 10.0,
4.0, 2.0, 3.0, 7.0, 25.0,
8.0, 6.0, 9.0, 4.0, 30.0,
9.0, 4.0, 3.0, 8.0, 45.0
};
プログラムは次のように書くことができます。
#define DIM 4
double a[DIM][DIM+1]={
2.0, 5.0, 8.0, 3.0, 10.0,
4.0, 2.0, 3.0, 7.0, 25.0,
8.0, 6.0, 9.0, 4.0, 30.0,
9.0, 4.0, 3.0, 8.0, 45.0
};
main()
{
int i,j,k;
double t;
/* 対角以外の要素を０にする */
for(i=0;i<DIM;i++){
for(j=0;j<DIM;j++){
if(j==i)continue;
t=a[j][i]/a[i][i];
for(k=0;k<=DIM;k++)a[j][k]-=a[i][k]*t;
}
disp(a); /* 動作が分かりやすいように経過を表示 */
getchar();
}
/* 対角の要素を１にする */
for(i=0;i<DIM;i++){
t=a[i][i];
for(j=0;j<=DIM;j++)a[i][j]/=t;
}
disp(a);
}
/* 行列の表示 */
disp(a)
double a[DIM][DIM+1];
{
int i,j;
for(i=0;i<DIM;i++){
for(j=0;j<=DIM;j++)printf("%10.4f ",a[i][j]);
printf("\n");
}
}
このプログラムを実行してみると、
2.0000 5.0000 8.0000 3.0000 10.0000
0.0000 -8.0000 -13.0000 1.0000 5.0000
0.0000 -14.0000 -23.0000 -8.0000 -10.0000
0.0000 -18.5000 -33.0000 -5.5000 0.0000
最初の列の対角以外の要素が０になっています。
2.0000 0.0000 -0.1250 3.6250 13.1250
0.0000 -8.0000 -13.0000 1.0000 5.0000
0.0000 0.0000 -0.2500 -9.7500 -18.7500
0.0000 0.0000 -2.9375 -7.8125 -11.5625
２番めの列の対角以外の要素が０になっています。
2.0000 0.0000 0.0000 8.5000 22.5000
0.0000 -8.0000 0.0000 508.0000 980.0000
0.0000 0.0000 -0.2500 -9.7500 -18.7500
0.0000 0.0000 0.0000 106.7500 208.7500
３番めの列の対角以外の要素が０になっています。
2.0000 0.0000 0.0000 -0.0000 5.8782
0.0000 -8.0000 0.0000 0.0000 -13.3958
0.0000 0.0000 -0.2500 0.0000 0.3162
0.0000 0.0000 0.0000 106.7500 208.7500
最後の列の対角以外の要素が０になっています。
1.0000 0.0000 0.0000 -0.0000 2.9391
0.0000 1.0000 0.0000 -0.0000 1.6745
0.0000 0.0000 1.0000 -0.0000 -1.2646
0.0000 0.0000 0.0000 1.0000 1.9555
対角の要素を１にしています。 ↑
これが解になる。
以外に簡単なプログラムで解くことができましたが、このアルゴリズムには問題があります。すべての行を総当たりで計算しますので、元の数が多くなると計算回数が飛躍的に増えてしまうのです。単純に計算すると、対角以外の要素を０にするループでは、ｎ元のときには ｎ＊ｎ＊（ｎ＋１） 回の計算を行ないますので、ｎが４のときには８０回、ｎが５のときには １５０回、ｎが１００のときには１０１万回の計算を行なわなければならない訳です。
続きはまた次回。
