連立１次方程式の解き方（掃出し法）

連立１次方程式

2x+ y+ +w=2

x+ 2z =0

x+ y+ z =1

y+2z+w=2

を解く。以下で、この方程式の係数と右辺の値を配列にして式変形を表している。この様にすれば計算機で解けるのであるが、原理を理解してプログラム化するために、紙の上で対応する式を書いて変形していくことを実行して下さい。

そのとき、紙の左半分に上の方程式、右半分に下の方程式(係数の部分は上と同じで、右辺の値が異なる)を書いて同じ式変形をすると一緒に解けることが見てとれます。

2x+ y+ +w=3

x+ 2z =-1

x+ y+ z =0

y+2z+w=1

(または、一般的に解く事を考えたら、右辺を定数を表す文字記号でやるとよい)

2x+ y+ +w=a

x+ 2z =b

x+ y+ z =c

y+2z+w=d

では、式変形を始めます。gap という数式処理システムでやっています。

gap>LogTo("lin_eq");

(この様に入力すると lin_eq というファイルにこの後の入出力が書込まれる)

gap> l:=[[2,1,0,1,2],[1,0,2,0,0],[1,1,1,0,1],[0,1,2,1,2]];

(この様に入力して左辺の記号 l を定める)

[ [ 2, 1, 0, 1, 2 ], [ 1, 0, 2, 0, 0 ], [ 1, 1, 1, 0, 1 ], [ 0, 1, 2, 1, 2 ] ]

(今、定めた l が表示される。方程式の係数と右辺を配列にしたものである)

gap> PrintArray(l);

[ [ 2, 1, 0, 1, 2 ],

[ 1, 0, 2, 0, 0 ],

[ 1, 1, 1, 0, 1 ],

[ 0, 1, 2, 1, 2 ] ]

(配列を見やすい形で表示して見る)

gap> l[1]:=1/l[1][1]*l[1];

[ 1, 1/2, 0, 1/2, 1 ]

(１式 l[1] の第１成分 l[1][1] の係数を １にするように１式を 1/l[1][1]倍する。

すると、今、定めた新しい l[1] が表示される)

gap> PrintArray(l);

[ [ 1, 1/2, 0, 1/2, 1 ],

[ 1, 0, 2, 0, 0 ],

[ 1, 1, 1, 0, 1 ],

[ 0, 1, 2, 1, 2 ] ]

gap> l[2]:=l[2]-l[2][1]*l[1];

[ 0, -1/2, 2, -1/2, -1 ]

(２式 l[2] から１式のl[2][1]倍を引いたものを新しい２式とする。これによって、

２式のxの係数 l[2][1] が0になる)

gap> PrintArray(l);

[ [ 1, 1/2, 0, 1/2, 1 ],

[ 0, -1/2, 2, -1/2, -1 ],

[ 1, 1, 1, 0, 1 ],

[ 0, 1, 2, 1, 2 ] ]

gap> l[3]:=l[3]-l[3][1]*l[1];

[ 0, 1/2, 1, -1/2, 0 ]

(３式 l[3] から１式のl[3][1]倍を引いたものを新しい３式とする)

gap> PrintArray(l);

[ [ 1, 1/2, 0, 1/2, 1 ],

[ 0, -1/2, 2, -1/2, -1 ],

[ 0, 1/2, 1, -1/2, 0 ],

[ 0, 1, 2, 1, 2 ] ]

(l[4][1]=0 であったので、これで第１列、つまり方程式のxの係数は１式以外は0に

なった)

gap> l[2]:=1/l[2][2]*l[2];

[ 0, 1, -4, 1, 2 ]

(２式 l[2] の第２成分 l[2][2] の係数を １にするように２式を 1/l[2][2]倍する)

gap> PrintArray(l);

[ [ 1, 1/2, 0, 1/2, 1 ],

[ 0, 1, -4, 1, 2 ],

[ 0, 1/2, 1, -1/2, 0 ],

[ 0, 1, 2, 1, 2 ] ]

gap> l[1]:=l[1]-l[1][2]*l[2];

[ 1, 0, 2, 0, 0 ]

(１式 l[1] から２式のl[1][2]倍を引いたものを新しい１式とする。これで１式のyの

係数 l[1][2] が0になった)

gap> PrintArray(l);

[ [ 1, 0, 2, 0, 0 ],

[ 0, 1, -4, 1, 2 ],

[ 0, 1/2, 1, -1/2, 0 ],

[ 0, 1, 2, 1, 2 ] ]

gap> l[3]:=l[3]-l[3][2]*l[2];

[ 0, 0, 3, -1, -1 ]

(３式 l[3] から２式のl[3][2]倍を引いたものを新しい３式とする。これで３式のyの

係数 l[3][2] が0になった)

gap> PrintArray(l);

[ [ 1, 0, 2, 0, 0 ],

[ 0, 1, -4, 1, 2 ],

[ 0, 0, 3, -1, -1 ],

[ 0, 1, 2, 1, 2 ] ]

gap> l[4]:=l[4]-l[4][2]*l[2];

[ 0, 0, 6, 0, 0 ]

(４式 l[4] から２式のl[4][2]倍を引いたものを新しい４式とする)

gap> PrintArray(l);

[ [ 1, 0, 2, 0, 0 ],

[ 0, 1, -4, 1, 2 ],

[ 0, 0, 3, -1, -1 ],

[ 0, 0, 6, 0, 0 ] ]

(第１列、つまり方程式のxの係数は１式以外は0のままで第２列、つまり方程式のyの

係数は２式以外は0になった)

gap> l[3]:=1/l[3][3]*l[3];

[ 0, 0, 1, -1/3, -1/3 ]

(３式 l[3] の第３成分 l[3][3] の係数を １にするように３式を 1/l[3][3]倍する)

gap> PrintArray(l);

[ [ 1, 0, 2, 0, 0 ],

[ 0, 1, -4, 1, 2 ],

[ 0, 0, 1, -1/3, -1/3 ],

[ 0, 0, 6, 0, 0 ] ]

gap> l[1]:=l[1]-l[1][3]*l[3];

[ 1, 0, 0, 2/3, 2/3 ]

(１式 l[1] から３式のl[1][3]倍を引いたものを新しい１式とする)

gap> PrintArray(l);

[ [ 1, 0, 0, 2/3, 2/3 ],

[ 0, 1, -4, 1, 2 ],

[ 0, 0, 1, -1/3, -1/3 ],

[ 0, 0, 6, 0, 0 ] ]

gap> l[2]:=l[2]-l[2][3]*l[3];

[ 0, 1, 0, -1/3, 2/3 ]

gap> PrintArray(l);

[ [ 1, 0, 0, 2/3, 2/3 ],

[ 0, 1, 0, -1/3, 2/3 ],

[ 0, 0, 1, -1/3, -1/3 ],

[ 0, 0, 6, 0, 0 ] ]

gap> l[4]:=l[4]-l[4][3]*l[3];

[ 0, 0, 0, 2, 2 ]

gap> PrintArray(l);

[ [ 1, 0, 0, 2/3, 2/3 ],

[ 0, 1, 0, -1/3, 2/3 ],

[ 0, 0, 1, -1/3, -1/3 ],

[ 0, 0, 0, 2, 2 ] ]

(これで、x、y、z の係数がそれぞれ１式、２式、３式以外は0になった)

gap> l[4]:=1/l[4][4]*l[4];

[ 0, 0, 0, 1, 1 ]

(最後に z についても４式を使って、１式、２式、３式から z を消去する)

gap> PrintArray(l);

[ [ 1, 0, 0, 2/3, 2/3 ],

[ 0, 1, 0, -1/3, 2/3 ],

[ 0, 0, 1, -1/3, -1/3 ],

[ 0, 0, 0, 1, 1 ] ]

gap> l[1]:=l[1]-l[1][4]*l[4];

[ 1, 0, 0, 0, 0 ]

gap> PrintArray(l);

[ [ 1, 0, 0, 0, 0 ],

[ 0, 1, 0, -1/3, 2/3 ],

[ 0, 0, 1, -1/3, -1/3 ],

[ 0, 0, 0, 1, 1 ] ]

gap> l[2]:=l[2]-l[2][4]*l[4];

[ 0, 1, 0, 0, 1 ]

gap> PrintArray(l);

[ [ 1, 0, 0, 0, 0 ],

[ 0, 1, 0, 0, 1 ],

[ 0, 0, 1, -1/3, -1/3 ],

[ 0, 0, 0, 1, 1 ] ]

gap> l[3]:=l[3]-l[3][4]*l[4];

[ 0, 0, 1, 0, 0 ]

gap> PrintArray(l);

[ [ 1, 0, 0, 0, 0 ],

[ 0, 1, 0, 0, 1 ],

[ 0, 0, 1, 0, 0 ],

[ 0, 0, 0, 1, 1 ] ]

これで方程式は

x =0

y =1

z =0

w =1

に変形されて解けた。

もし、与えられた方程式の１式の最初の係数が0ならば、そこが0にならない様に方程式の式の並べ方を変える。次の l2 はその様な方程式の例の配列です。

gap> l2:=[[0,1,2,1,2],[1,0,2,0,0],[1,1,1,0,1],[2,1,0,1,2]];

[ [ 0, 1, 2, 1, 2 ], [ 1, 0, 2, 0, 0 ], [ 1, 1, 1, 0, 1 ], [ 2, 1, 0, 1, 2 ] ]

gap> PrintArray(l2);

[ [ 0, 1, 2, 1, 2 ],

[ 1, 0, 2, 0, 0 ],

[ 1, 1, 1, 0, 1 ],

[ 2, 1, 0, 1, 2 ] ]

gap> l2:=l2{[4,2,3,1]};

[ [ 2, 1, 0, 1, 2 ], [ 1, 0, 2, 0, 0 ], [ 1, 1, 1, 0, 1 ], [ 0, 1, 2, 1, 2 ] ]

(１式と４式を取り換えた。前の式変形で、1/l[2][2] を計算するときに l[2][2] が0 ならば、３式または４式のyの係数 l[3][2], l[4][2] の0でない式と２式を交換する。 1/l[3][3] を計算するときに l[3][3] が0のときは３式と４式を取り換える。この様 にしても0がさけられないときは、実は、解は唯１個にはならないのである)

gap> PrintArray(l2);

[ [ 2, 1, 0, 1, 2 ],

[ 1, 0, 2, 0, 0 ],

[ 1, 1, 1, 0, 1 ],

[ 0, 1, 2, 1, 2 ] ]

gap> LogTo();

(これで入出力の記録が終りになります)

いろいろな方程式で上の操作を実行して見て下さい。1/l[n][n] をするときに障害が

でると思います。次の例は大変極端な場合です。

w = 3

x = -1

y = 2

z = 1

この場合は １式と２式の交換 l{[2,1,3,4]} をして、さらに２式と３式の交換l{[1,3,2,4]} をして、さらに３式と４式の交換 l{[1,2,4,3]} をすることになります。人が見れば始めから解く必要のない方程式であることは明らかなのですが、プログラムを作って解くときには、、、、何かいい知恵はありますか？