LUa_ecomposition.zip三角関数、逆三角関数、その他関数、 連立一次方程式の解法 に戻る
連立一次方程式の解法６
連立一次方程式の解法の続きです。
今回のものは、C言語辞典にあるLU分解をDelphi用に変換し、判定を追加したものです。
基本部分は、連立方程式の解法３と同じですが、weightの計算とpivot操作があり、行列式の値が求められるようになっています。
行列式の値から、演算結果が正しいかどうかの判別を追加しようと思いましたが、値の設定がおもわしくなく、取りやめました。
プログラム
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
darray = array of array of extended;
sarray = array of extended;
iarray = array of integer;
type
TForm1 = class(TForm)
BitBtn1: TBitBtn;
Memo1: TMemo;
StringGrid1: TStringGrid;
StringGrid2: TStringGrid;
LabeledEdit1: TLabeledEdit;
BitBtn2: TBitBtn;
LabeledEdit2: TLabeledEdit;
CheckBox1: TCheckBox;
RadioGroup1: TRadioGroup;
procedure BitBtn1Click(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure BitBtn2Click(Sender: TObject);
procedure RadioGroup1Click(Sender: TObject);
private
{ Private 宣言 }
procedure GridSet(Sin, Colin: integer); // グリッドサイズ変更
function gauss(n: integer; var a: darray; var b, x: sarray): extended;
function lu(n: integer; var a: darray; var ip: iarray): extended;
procedure solve(n : integer; var a: darray; var b, x: sarray; var ip: iarray);
procedure testset;
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, 34, 24, 32),
( 2, 7, 0, 0),
(24, 17, 0, 0),
( 7, 0, 0, 0));
mattestD : 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));
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, Asave : darray; // A行列配列
b, x : sarray;
epsilon : extended;
// グリッドの設定
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;
// lu分解
// n 行列数
// a 行列の値
// ip 行順
// 戻り値 行列の値
function TForm1.lu(n: integer; var a: darray; var ip: iarray): extended;
var
i, j, k, ii, ik: integer;
t, u : extended;
weight : sarray;
begin
setlength(weight, n);
result := 1; // <> 0 成功 0 失敗
j := 0;
for k := 0 to n - 1 do begin // 各行 重み計算
u := 0; // その行の絶対値最大の要素を求める
for j := 0 to n - 1 do begin
t := abs(a[k, j]);
if t > u then u := t; // 最大値
end;
j := 0;
if u = 0 then exit; // 最大値が0 なら行列はLU分解できない 同じ列が全てゼロ
weight[k] := 1 / u; // 最大値の逆数
end;
for k := 0 to n - 1 do begin //
u := -1;
for i := k to n - 1 do begin // kからの行
ii := ip[i]; // 重み×絶対値 が最大の行を見つける
t := abs(a[ii, k]) * weight[ii];
if t > u then begin
u := t;
j := i; // j 最大の行
end;
end;
ik := ip[j];
if j <> k then begin // 同じ行でなかったら行をを交換
ip[j] := ip[k];
ip[k] := ik;
result := -result;
end;
u := a[ik, k]; // 対角成分 行が入れ替わっているのでik<>kあり
if u = 0 then begin
result := 0;
exit; // 対角にゼロが有ったら計算出来ないので終了
end;
result := result * u;
for i := k + 1 to n - 1 do begin // ガウスの消去法
ii := ip[i];
a[ii, k] := a[ii, k] / u; // L
t := a[ii, k];
for j := k + 1 to n - 1 do a[ii, j] := a[ii, j] - t * a[ik, j]; // U
j := 0;
end;
end;
end;
// LUを使用して連立方程式を解く
procedure TForm1.solve(n : integer; var a: darray; var b, x: sarray; var ip: iarray);
var
i, j, ii: integer;
t : extended;
begin
for i := 0 to n - 1 do begin // ガウスの消去法の残り 前進代入
ii := ip[i];
t := b[ii];
for j := 0 to i - 1 do t := t - a[ii, j] * x[j];
x[i] := t;
end;
for i := n - 1 downto 0 do begin // 後退代入
t := x[i];
ii := ip[i];
for j := i + 1 to n - 1 do t := t - a[ii, j] * x[j];
x[i] := t / a[ii, i];
end;
end;
// ガウス法
function TForm1.gauss(n: integer; var a: darray; var b, x: sarray): extended;
var
ip : iarray;
det : extended;
i : integer;
LL : string;
begin
setlength(ip, n);
for i := 0 to n - 1 do ip[i] := i; // 行No初期化
det := lu(n, a, ip); // LU分解
if det = 0 then
MessageDlg('対角にゼロが発生し計算出来ません。', mtInformation,[mbOk], 0);
if det <> 0 then begin // LU分解出来たら
solve(n, a, b, x, ip); // LUを使用して連立方程式を解く
LL := '計算順 ';
for i := 0 to n - 2 do begin
LL := LL + inttostr(ip[i] + 1) + ', ';
end;
LL := LL + inttostr(ip[n -1] + 1);
memo1.Lines.Append(LL);
end;
result := det; // 戻り値 <> 0 成功 0 失敗
end;
// 計算と検算
procedure TForm1.BitBtn1Click(Sender: TObject);
var
i, j, ch: integer;
det, aer : extended;
y : sarray;
allowable_errorr : extended;
maxb, minb, bc : extended;
acmax : sarray;
acmin : sarray;
begin
setlength(matA, Nl, Nl);
setlength(Asave, Nl, Nl);
setlength(b, Nl);
setlength(x, Nl);
setlength(y, Nl);
setlength(acmax, Nl);
setlength(acmin, 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;
for i := 0 to Nl - 1 do // matAを Asaveに保存
for j := 0 to Nl - 1 do Asave[i, j] := matA[i, j];
// 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;
// 検算許容誤差
val(labelededit2.Text, allowable_errorr, ch);
if ch <> 0 then begin
MessageDlg('検算誤差許容値の入力値の' + intTostr(ch) + '文字目に誤りが有ります。', mtInformation,[mbOk], 0);
exit;
end;
// A行列の最大値最小値検索
for i := 0 to Nl -1 do begin
acmax[i] := 1;
acmin[i] := 1;
if matA[i, 0] <> 0 then acmax[i] := abs(matA[i, 0]);
if matA[i, 0] <> 0 then acmin[i] := abs(matA[i, 0]);
for j := 1 to Nl - 1 do begin
if (abs(matA[i, j]) > acmax[i]) and (matA[i, j] <> 0) then acmax[i] := abs(matA[i, j]);
if (abs(matA[i, j]) < acmin[i]) and (matA[i, j] <> 0) then acmin[i] := abs(matA[i, j]);
end;
end;
// B行列とA行列の比最小値と最大値検索
maxb := abs(b[0]) / acmax[0];
minb := abs(b[0]) / acmin[0];
for i := 1 to Nl - 1 do begin
if abs(b[i]) / acmax[i] > maxb then maxb := abs(b[i]) / acmax[i];
if abs(b[i]) / acmin[i] < minb then minb := abs(b[i]) / acmin[i];
end;
// bの値の範囲が演算桁数の範囲を越えているか判定
bc := minb;
if (minb = 0) and (maxb > 0) then
if maxb >= 1 then bc := 1 / maxb
else bc := maxb;
if (minb > 0) and (maxb > 0) then bc := minb / maxb;
if (minb = 0) and (maxb = 0) then bc := 1;
if checkbox1.Checked = false then
if bc < epsilon then begin
MessageDlg('bの入力値の値が演算の範囲を超えています。', mtInformation,[mbOk], 0);
memo1.Clear;
exit;
end;
// 検算結果の判定値設定
allowable_errorr := allowable_errorr * sqrt(bc);
memo1.Clear;
det := gauss(Nl, matA, b, x); // Gausas で Axn = b で解く det 行列式の値
if det = 0 then exit;
memo1.Lines.Append('行列式の値 ' + floatTostr(det));
// 検算差分の表示
memo1.Lines.Append('行No 差分');
for i := 0 to Nl - 1 do begin
aer := 0;
for j := 0 to Nl - 1 do begin
aer := aer + Asave[i, j] * x[j];
end;
y[i] := b[i] - aer;
memo1.Lines.Append('No' + intTostr(i + 1) + ' ' + floatTostr(y[i]));
end;
j := 0;
for i := 0 to Nl - 1 do begin
if abs(b[i]) > allowable_errorr then aer := y[i] / b[i]
else aer := y[i];
if abs(aer) > allowable_errorr then inc(j);
end;
if checkbox1.Checked = false then begin
if j = 0 then memo1.Lines.Append('判定 OK')
else memo1.Lines.Append('判定 NG');
if j > 0 then exit;
end;
for i := 0 to Nl - 1 do memo1.Lines.Append('x' + intTostr(i + 1) + '= ' + floatTostr(x[i]));
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.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.RadioGroup1Click(Sender: TObject);
begin
testset;
end;
// 初期値の設定
procedure TForm1.FormCreate(Sender: TObject);
var
a : extended;
begin
testset;
epsilon := 1;
repeat
epsilon := epsilon / 2;
a := 1 - epsilon;
until a = 1;
epsilon := epsilon * 2;
end;
end.
LUa_ecomposition.zip
三角関数、逆三角関数、その他関数、 連立一次方程式の解法 に戻る