LU_test.zip三角関数、逆三角関数、その他関数、 連立一次方程式の解法 に戻る
連立一次方程式の解法２
連立一次方程式の解法１の続きです。
対角にゼロが入るかどうかをLU分解に先だって三角行列を作成して、行列の値から判別します。
対角にゼロが入らなければ、ＬＵ分解を行いますが、LU分解が出来ない場合は、LU分解が出来るまで、行の変更を行います。
今回の場合は、行の変更を、行の組み合わせ全てについて先に作成しておいて、順次組み合わせを変更して、解が求まるまで繰り返します。
解が求まらない場合は、全ての組み合わせの計算をするので解の無い方程式の行列式となります。
組み合わせの数は、行の数をｎとすると、ｎの階乗となります。
５×５の場合で、２４、６×６のの場合１２０もの組み合わせとなるのですが、PCでの実行なので、さほど時間は掛かりません。
三角行列の行列式の値から対角にゼロが入るかどうか判別していますが、三角行列の作成時除算が入るので、 演算誤差によ数学的には行列式の値がゼロになる値が演算上必ずゼロになる保証は有りません。
LU分解時も同じなので、最終的には、検算による判別が必要ですが、此処のプログラムでは、最終判別を省略しています。
プログラム
unit Main;
interface
uses
Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, Vcl.Buttons, Vcl.Grids, system.UITypes,
Vcl.ExtCtrls;
type
TForm1 = class(TForm)
BitBtn1: TBitBtn;
Memo1: TMemo;
StringGrid1: TStringGrid;
StringGrid2: TStringGrid;
LabeledEdit1: TLabeledEdit;
BitBtn2: TBitBtn;
RadioGroup1: TRadioGroup;
procedure BitBtn1Click(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure BitBtn2Click(Sender: TObject);
procedure RadioGroup1Click(Sender: TObject);
private
{ Private 宣言 }
function LUecomposition: byte;
procedure GridSet(Sin, Colin: integer); // グリッドサイズ変更
procedure rowmatinc(k: integer);
procedure testset;
function triangular_matrix(var z : integer): extended;
procedure permutation;
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
var
Nltest : integer = 4; // テストデーター行列数
// テストデーター
mattestA : array[0..3] of array[0..3] of extended = ((8, 16, 24, 32),
(2, 7, 12, 17),
(6, 17, 32, 59),
(7, 22, 46,105));
mattestB : array[0..3] of array[0..3] of extended = ((8, 16, 24, 32),
(2, 7, 12, 17),
(2, 17, 32, 0),
(7, 0, 0, 105));
mattestC : array[0..3] of array[0..3] of extended = ((8, 16, 24, 0),
(2, 7, 12, 0),
(2, 17, 32, 0),
(7, 0, 0, 105));
mattestD : array[0..3] of array[0..3] of extended = ((0, 16, 24, 32),
(2, 7, 12, 17),
(6, 17, 32, 0),
(7, 0, 0, 105));
mattestE : array[0..3] of array[0..3] of extended = ((0, 0, 0, 32),
(0, 0, 12, 17),
(6, 17, 32, 0),
(7, 0, 0, 105));
btestA : array[0..3] of extended = (160, 70, 198, 291);
btestB : array[0..3] of extended = (150, 60, 10, 40);
Nl : integer; // 行列数
matA : array of array of extended; // A行列配列
matL : array of array of extended; // L行列配列
matU : array of array of extended; // U行列配列
lu : array of array of extended;
b : array of extended;
ba : array of extended;
a : array of array of extended;
perNo : array of array of integer; // 全行No配列
rowNo : array of integer;
exchange : integer;
diagonal : integer;
// グリッドの設定
procedure TForm1.GridSet(Sin, Colin: integer); // グリッドサイズ変更
var // データー行の数は10迄
H, W, i :integer; // それ以上はスクロールする
begin
inc(Sin);
inc(Colin);
StringGrid1.ColCount := Colin;
StringGrid1.RowCount := Sin;
if Sin <= 11 then // ストリンググリッドの大きさ設定
begin // 固定行を含め11行又はそれ以下の場合
H := (StringGrid1.DefaultRowHeight + 1) * Sin;
StringGrid2.ScrollBars := ssNone; // スクロールバー表示
end
else
begin // 固定行を含め11行より多い場合
H := (StringGrid1.DefaultRowHeight + 1) * (10 + 1);
StringGrid2.ScrollBars := ssVertical; // スクロールバー表示
end;
if Colin <= 11 then // ストリンググリッドの大きさ設定
begin // 固定行を含め11行又はそれ以下の場合
W := (StringGrid1.DefaultColWidth + 1) * Colin;
StringGrid1.ScrollBars := ssNone; // スクロールバーなし
end
else
begin // 固定行を含め11行より多い場合
W := (StringGrid1.DefaultColWidth + 1) * (10 + 1);
StringGrid1.ScrollBars := ssBoth; // スクロールバー表示
end;
StringGrid1.ClientWidth := W;
StringGrid1.ClientHeight := H;
for i := 1 to Sin - 1 do StringGrid1.Cells[0,i] := 'a' + intTostr(i) + '*';
for i := 1 to Sin - 1 do StringGrid1.Cells[i,0] := 'a*' + intTostr(i);
StringGrid2.RowCount := Sin;
StringGrid2.ClientHeight := H + 1;
StringGrid2.ClientWidth := (StringGrid2.DefaultColWidth + 1) * 2;
for i := 1 to Sin - 1 do StringGrid2.Cells[0, i] := intTostr(i);
StringGrid2.Cells[1, 0] := 'b*=';
end;
// 行の入れ替え 対角にゼロが発生した時に行の入れ替えを行います。
procedure TForm1.rowmatinc(k: integer);
var
i : integer;
begin
for i := 0 to Nl - 1 do rowNo[i] := perNo[k, i];
end;
// 行列式の値
// 単に対角にゼロがあるかどうかの確認
// 対角の計算値がゼロになのに演算誤差でゼロにならないときあり
function TForm1.triangular_matrix(var z : integer): extended;
var
i, j, k, m, q: integer;
buf, det: extended;
c : array of array of extended;
p : array of integer;
begin
setlength(c, Nl, Nl);
setlength(p, Nl);
for i := 0 to Nl - 1 do begin
p[i] := i;
for j := 0 to Nl - 1 do
c[i, j] := matA[i, j];
end;
z := 0;
det := 1;
// 三角行列の作成
for i := 0 to Nl -1 do
for j := 0 to Nl -1 do
if i < J then begin
if c[i, i] = 0 then begin
for k := i + 1 to Nl - 1 do
if c[k, i] <> 0 then begin
det := -det;
for m := 0 to Nl - 1 do begin
buf := c[i, m];
c[i, m] := c[k, m];
c[k, m] := buf;
end;
q := p[i];
p[i] := p[k];
p[k] := q;
end;
end
else begin
buf := c[j, i] / c[i, i];
for k := 0 to Nl - 1 do c[j, k] := c[j, k] - c[i, k] * buf;
end;
end;
// 対角行列の積
for i := 0 to Nl - 1 do begin
det := det * c[i, i];
if c[i, i] = 0 then z := p[i] + 1;
end;
result := det;
memo1.Lines.Append('行列式の値 ' + floattostr(det));
end;
// 全行No組み合わせ配列作成
procedure TForm1.permutation;
var
v : array of integer;
i, n : integer;
procedure swap(var pa, pb : integer);
var
tmp : integer;
begin
tmp := pa;
pa := pb;
pb := tmp;
end;
procedure reverse(first, last: integer);
var
lam : integer;
begin
lam := last;
dec(lam);
while (first <> last) and (first <> lam) do begin
swap(v[first], v[lam]);
inc(first);
last := lam;
dec(lam);
end;
end;
function next_permutation(first, last: integer): boolean;
var
i, j, k : integer;
begin
if first = last then begin
result := false;
exit;
end;
if first + 1 = last then begin
result := false;
exit;
end;
i := last - 1;
while true do begin
j := i;
dec(i);
if v[i] < v[j] then begin
k := last;
dec(k);
while not(v[i] < v[k]) do dec(k);
swap(v[i], v[k]);
reverse(j, last);
result := true;
exit;
end;
if i = first then begin
reverse(j, last);
result := false;
exit;
end;
end;
end;
begin
setlength(v, Nl);
for i := 0 to Nl -1 do v[i] := i;
n := 0;
repeat
for i := 0 to Nl - 1 do perNo[n, i] := v[i];
inc(n);
until not next_permutation(0, Nl);
end;
// LU分解による連立方程式の解法
function TForm1.LUecomposition: byte;
const
sp = ' ';
var
i, j, k, n : integer;
sum : extended;
aa : array of array of extended;
l, u, x, y : array of extended;
LL, spl : string;
begin
setlength(a, Nl, Nl);
setlength(aa, Nl, Nl);
setlength(l, Nl);
setlength(u, Nl);
setlength(x, Nl);
setlength(y, Nl);
// matデーターを編集用にコピー matA, B
for j := 0 to Nl - 1 do begin
ba[j ] := b[rowNo[j]];
for k := 0 to Nl - 1 do a[j, k] := matA[rowNo[j], k];
end;
// 編集用をバックアップににコピー
for j := 0 to Nl - 1 do
for k := 0 to Nl - 1 do
aa[j, k] := a[j, k];
// matUの対角を1に設定
for i := 0 to Nl - 1 do
for j := 0 to Nl - 1 do
if i = j then begin
matU[i, j] := 1;
matL[i, j] := 0;
end
else begin
matU[i, j] := 0;
matL[i, j] := 0;
end;
// LU分解
for i := 0 to Nl -1 do begin
n := Nl - i - 1;
// matA先頭の値をmatLの対角にセット
// 対角がゼロだったら終了
if a[0, 0] = 0 then begin
break;
end;
matL[i, i] := a[0, 0];
// MatLの作成
for j := 0 to n - 1 do begin
l[j] := a[j + 1, 0];
matL[j + i + 1, i] := l[j];
end;
// matUの作成
for j := 0 to n - 1 do begin
u[j] := a[0, j + 1] / a[0, 0];
matU[i, j + i + 1] := u[j];
end;
// luを求める
for j := 0 to n - 1 do begin
for k := 0 to n - 1 do begin
lu[j, k] := l[j] * u[k];
end;
end;
// aを求める
for j := 0 to n - 1 do begin
for k := 0 to n - 1 do begin
a[j, k] := a[j + 1, k + 1] - lu[j, k];
end;
end;
end;
diagonal := Nl;
// 対角にゼロが発生した場合演算を中止し'
// 行の入れ替えを試みます
if a[0, 0] = 0 then begin
result := 1;
diagonal := i;
exit;
end;
// LU行列で連立方程式を解きます
for i := 0 to Nl - 1 do begin
sum := 0;
for k := 0 to i - 1 do sum := sum + matL[i, k] * y[k];
y[i] := (ba[i] - sum) / matL[i, i];
end;
// 後退代入
for i := Nl - 1 downto 0 do begin
sum := 0;
for k := i + 1 to Nl - 1 do sum := sum + matU[i, k] * x[k];
x[i] := y[i] - sum;
end;
// 各行代入答えの差分計算
for j := 0 to Nl - 1 do begin
sum := 0;
for i := 0 to Nl - 1 do begin
sum := sum + aa[j, i] * x[i];
end;
y[j] := ba[j] - sum;
end;
// 演算行No表示
LL := '計算順 ';
for i := 0 to Nl - 1 do begin
if i < NL -1 then LL := LL + intTostr(rowNo[i] + 1) + ','
else LL := LL + intTostr(rowNo[i] + 1);
end;
memo1.Lines.Append(LL);
// 行番号表示
memo1.Lines.Append('行 No 検算による差分');
for k := 0 to Nl - 1 do begin
LL := 'No = ' + inttostr(k + 1);
n := 25 - length(LL);
spl := '';
for j := 1 to n do spl := spl + sp;
memo1.Lines.Append(LL + spl + floattostr(y[rowNo[k]]));
end;
// xの値表示
for k := 0 to Nl - 1 do begin
memo1.Lines.Append('X' + inttostr(k + 1) + '=' + floatTostr(x[k]));
end;
result := 0;
end;
procedure TForm1.BitBtn1Click(Sender: TObject);
var
i, j, ch, k : integer;
LL : string;
begin
setlength(matA, Nl, Nl);
setlength(matL, Nl, Nl);
setlength(b, Nl);
setlength(ba, Nl);
setlength(lu , Nl, Nl);
setlength(matU, Nl, Nl);
setlength(rowNo, Nl);
// matAにグリッド1から値読み込み
for i := 0 to Nl - 1 do
for j := 0 to Nl - 1 do begin
val(StringGrid1.Cells[i + 1,j + 1],matA[j,i], ch);
if ch <> 0 then begin
MessageDlg('a' + intTostr(j + 1) + ',' + intTostr(i + 1) + ' の入力値の' + intTostr(ch) + '文字目に誤りが有ります。', mtInformation,[mbOk], 0);
exit;
end;
end;
// bにグリッド2から値読み込み
for i := 0 to Nl - 1 do begin
val(StringGrid2.Cells[1, i + 1], b[i], ch);
if ch <> 0 then begin
MessageDlg('b' + intTostr(i + 1) + 'の入力値の' + intTostr(ch) + '文字目に誤りが有ります。', mtInformation,[mbOk], 0);
exit;
end;
end;
memo1.Clear;
// 対角にゼロが有ったら終了
if triangular_matrix(k) = 0 then begin
MessageDlg('配列Aのaにゼロがあるか'
+ #13#10 + 'LU分解時ゼロが発生します。', mtInformation,[mbOk], 0);
exit;
end;
k := 2;
for i := 3 to NL do k := k * i; // 行組み合わせ総数
setlength(perNo, k, Nl); // 全行組み合わせ配列確保
permutation; // 全行組み合わせ作成
for i := 0 to Nl - 1 do rowNo[i] := i; // 計算準初期化
ch := k; // 最大ループ数
// Lu分解連立方程式の解法 解法できない場合あり
exchange := 0;
j := 0; // 解法できた場合 0 出来ない場合 1 対角値0の場合 2
// Lu分解が正常に終了するまで繰り返す
// k = 0 は初期設定値なので1～
for k := 1 to ch do begin
j := LUecomposition; // LU分解
if j = 1 then begin // 対角にゼロが発生
rowmatinc(k); // 行の移動
inc(exchange);
end;
if j = 0 then break; // 解法で終了
end;
// 解法出来なかったら
if diagonal < Nl then
MessageDlg('LU分解中に対角ゼロが発生しました、入力値を見直して下さい。', mtInformation,[mbOk], 0);
memo1.Lines.Append('再計算 Row= ' + intTostr(exchange));
if j = 1 then begin
LL := '計算順 ';
for i := 0 to Nl - 1 do begin
if i < NL -1 then LL := LL + intTostr(rowNo[i] + 1) + ','
else LL := LL + intTostr(rowNo[i] + 1);
end;
memo1.Lines.Append(LL);
end;
end;
// テストデーターの選択設定
procedure TForm1.testset;
var
i, j : integer;
begin
GridSet(Nltest, Nltest);
Nl := Nltest;
for i := 1 to Nltest do
for j := 1 to Nltest do begin
case radiogroup1.ItemIndex of
0: StringGrid1.Cells[i,j] := floatTostr(mattestA[j-1, i-1]);
1: StringGrid1.Cells[i,j] := floatTostr(mattestB[j-1, i-1]);
2: StringGrid1.Cells[i,j] := floatTostr(mattestC[j-1, i-1]);
3: StringGrid1.Cells[i,j] := floatTostr(mattestD[j-1, i-1]);
4: StringGrid1.Cells[i,j] := floatTostr(mattestE[j-1, i-1]);
end;
end;
for i := 1 to Nltest do begin
case radiogroup1.ItemIndex of
0: StringGrid2.Cells[1,i] := floatTostr(btestA[i-1]);
1, 2, 3, 4 : StringGrid2.Cells[1,i] := floatTostr(btestB[i-1]);
end;
end;
LabeledEdit1.Text := intTostr(Nltest);
memo1.Clear;
end;
// グリッドサイズの変更
procedure TForm1.BitBtn2Click(Sender: TObject);
var
ch : integer;
begin
val(LabeledEdit1.Text, Nl, ch);
if ch <> 0 then begin
MessageDlg('N数の入力値の' + intTostr(ch) + '文字目に誤りが有ります。', mtInformation,[mbOk], 0);
exit;
end;
GridSet(Nl, Nl);
end;
// テストデーターセット
procedure TForm1.RadioGroup1Click(Sender: TObject);
begin
testset;
end;
// 初期値の設定
procedure TForm1.FormCreate(Sender: TObject);
begin
testset;
end;
end.
LU_test.zip
三角関数、逆三角関数、その他関数、 連立一次方程式の解法 に戻る