LU_Bigdecimal.zip三角関数、逆三角関数、その他関数、 連立一次方程式の解法 に戻る
連立一次方程式の解法７
連立一次方程式の解法の続きです。
今回は、連立一次方程式の解法５をBigDecimalに変更したものです。
行列式の計算の中には、特に問題になる演算はなく、ExtendedをBigDecimalに変更すれば正しい演算結果を得られます。
単に、計算結果の有効桁数が増えるだけです。
有効桁数の入力を60迄に制限していますが、単に表示の問題だけで制限しているので変更しても問題ありません。
プログラム
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, Velthuis.BigIntegers, Velthuis.BigDecimals;
type
TForm1 = class(TForm)
BitBtn1: TBitBtn;
Memo1: TMemo;
StringGrid1: TStringGrid;
StringGrid2: TStringGrid;
LabeledEdit1: TLabeledEdit;
BitBtn2: TBitBtn;
RadioGroup1: TRadioGroup;
LabeledEdit2: TLabeledEdit;
LabeledEdit3: TLabeledEdit;
CheckBox1: TCheckBox;
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 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 = ((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; // 行列数
aa : array of array of bigdecimal; // 演算用a配列
a : array of array of bigdecimal; // 演算用a配列
b : array of bigdecimal; // 演算用b行列
p : array of integer; // 行移動
dpcs : integer;
allowable_errorr : bigdecimal;
// グリッドの設定
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分解
function TForm1.LUecomposition: byte;
const
sp = ' ';
var
i, j, k, n, pivot, pt: integer;
sum, tmp, amax : bigdecimal;
x, y : array of bigdecimal;
LL, spl : string;
begin
setlength(a, Nl, Nl);
setlength(x, Nl);
setlength(y, Nl);
setlength(p, Nl);
memo1.Clear;
// matA, bデーターを編集用bigdecimalにコピー
// aの値は変化してしまうため
for j := 0 to Nl - 1 do begin
for k := 0 to Nl - 1 do begin
a[j, k] := aa[j, k];
end;
end;
// 行入れ替えno初期化
for i := 0 to Nl - 1 do p[i] := i; // 行番号0～N - 1
// LU分解 前進消去
for k := 0 to Nl - 2 do begin
pivot := k; // ピボットの選択
amax := bigdecimal.abs(a[pivot, pivot]); // K列の行の最大値検索
for i := pivot + 1 to Nl - 1 do begin // ピボットの次の行から
if bigdecimal.abs(a[i, k]) > amax then begin
pivot := i;
amax := bigdecimal.abs(a[pivot, k]);
end;
end;
if pivot <> k then begin // 最大値行が違ったら入れ替え
for i := 0 to Nl - 1 do begin
tmp := a[k, i];
a[k, i] := a[pivot, i];
a[pivot, i] := tmp;
end;
pt := p[k]; // 行番号交換
p[k] := p[pivot];
p[pivot] := pt;
end;
// LU decomposition right-looking algorithm
for i := k + 1 to Nl - 1 do begin
if a[k, k] = 0 then begin // 対角にゼロが発生したら終了
result := 1;
exit;
end;
a[i, k] := a[i, k] / a[k, k]; // L
for j := k + 1 to Nl - 1 do a[i, j] := a[i, j] - a[i, k] * a[k, j]; // U
end;
end;
// bの値セット
for i := 0 to Nl - 1 do x[i] := b[p[i]]; // 行入れ替えnoによりxにbの値セット
// 前進代入
for i := 1 to Nl - 1 do
for j := 0 to i - 1 do x[i] := x[i] - a[i, j] * x[j];
// 後退代入
for i := Nl - 1 downto 0 do begin
for j := i + 1 to Nl - 1 do x[i] := x[i] - a[i, j] * x[j];
if a[i, i] = 0 then begin // 対角ゼロが発生したら終了
result := 1;
exit;
end;
x[i] := x[i] / a[i, i];
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]; // 元のaa[]の値を使用します。
end;
y[j] := b[j] - sum; // 差分
end;
LL := '計算順 ';
for i := 0 to Nl - 1 do begin
if i < NL - 1 then LL := LL + intTostr(p[i] + 1) + ','
else LL := LL + intTostr(p[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 := 10 - length(LL);
spl := '';
for j := 1 to n do spl := spl + sp;
memo1.Lines.Append(LL + spl + y[k].ToString);
end;
// 誤差の確認 誤差が大きい場合演算途中対角の値がゼロに近い値になっています。
j := 0; // 誤差大の数
for k := 0 to Nl - 1 do begin
if bigdecimal.abs(b[k]) > allowable_errorr then amax := y[k] / b[k] // 誤差比率
else amax := y[k];
if bigdecimal.abs(amax) > allowable_errorr then inc(j);
end;
if checkbox1.Checked = false then
if j > 0 then begin
memo1.Lines.Append('検算結果 NG');
result := 1;
exit;
end
else
memo1.Lines.Append('検算結果 OK');
// xの値表示
for k := 0 to Nl - 1 do begin
x[k] := x[k].RemoveTrailingZeros(0);
memo1.Lines.Append('X' + inttostr(k + 1) + '=' + x[k].ToString);
end;
result := 0;
end;
procedure TForm1.BitBtn1Click(Sender: TObject);
var
i, j, ch: integer;
matA : array of array of extended; // A行列配列 入力チェック用 演算には使用されない
matB : array of extended; // A行列配列 入力チェック用 演算には使用されない
maxb, minb, bc, epsilon : bigdecimal;
LL : string;
acmax : array of bigdecimal;
acmin : array of bigdecimal;
begin
setlength(matA, Nl, Nl);
setlength(aa, Nl, Nl);
setlength(matB, Nl);
setlength(b, 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;
// bにグリッド2から値読み込み
for i := 0 to Nl - 1 do begin
val(StringGrid2.Cells[1, i + 1], matB[i], ch);
if ch <> 0 then begin
MessageDlg('b' + intTostr(i + 1) + 'の入力値の' + intTostr(ch) + '文字目に誤りが有ります。', mtInformation,[mbOk], 0);
exit;
end;
end;
// bigdecimalの配列に読み込み
for i := 0 to Nl - 1 do begin
b[i] := StringGrid2.Cells[1, i + 1];
for j := 0 to Nl - 1 do begin
aa[j, i] := StringGrid1.Cells[i + 1,j + 1];
end;
end;
val(labelededit2.Text, dpcs, ch);
if ch <> 0 then begin
MessageDlg('有効桁数入力値' + intTostr(ch) + '文字目に誤りが有ります。', mtInformation,[mbOk], 0);
exit;
end;
if dpcs < 5 then begin
MessageDlg('有効桁数の値が小さすぎます。', mtInformation,[mbOk], 0);
exit;
end;
if dpcs > 60 then begin
MessageDlg('有効桁数のは60迄にしてください。', mtInformation,[mbOk], 0);
exit;
end;
// 誤差の許容範囲読み込み
val(labelededit3.Text, matB[0], ch);
if ch <> 0 then begin
MessageDlg('許容誤差の' + intTostr(ch) + '文字目に誤りが有ります。', mtInformation,[mbOk], 0);
exit;
end;
if matB[0] > 10 then begin
MessageDlg('許容誤差の値が大きすぎます。', mtInformation,[mbOk], 0);
exit;
end;
if matB[0] <= 0 then begin
MessageDlg('許容誤差はゼロより大きくして下さい。', mtInformation,[mbOk], 0);
exit;
end;
allowable_errorr := labelededit3.Text;
// epsilonを有効桁数から設定
LL := '1e-' + intTostr(dpcs - 1);
epsilon := LL;
// A行列の最大値最小値検索
for i := 0 to Nl -1 do begin
acmax[i] := 1;
acmin[i] := 1;
if aa[i, 0] <> 0 then acmax[i] := bigdecimal.abs(aa[i, 0]);
if aa[i, 0] <> 0 then acmin[i] := bigdecimal.abs(aa[i, 0]);
for j := 1 to Nl - 1 do begin
if (bigdecimal.abs(aa[i, j]) > acmax[i]) and (aa[i, j] <> 0) then acmax[i] := bigdecimal.abs(aa[i, j]);
if (bigdecimal.abs(aa[i, j]) < acmin[i]) and (aa[i, j] <> 0) then acmin[i] := bigdecimal.abs(aa[i, j]);
end;
end;
// B行列とA行列の比最小値と最大値検索
maxb := bigdecimal.abs(b[0]) / acmax[0];
minb := bigdecimal.abs(b[0]) / acmin[0];
for i := 1 to Nl - 1 do begin
if bigdecimal.abs(b[i]) / acmax[i] > maxb then maxb := bigdecimal.abs(b[i]) / acmax[i];
if bigdecimal.abs(b[i]) / acmin[i] < minb then minb := bigdecimal.abs(b[i]) / acmin[i];
end;
if minb > 1 then minb := 1 / minb;
if maxb < 1 then maxb := 1 / maxb;
// 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 * bigdecimal.Sqrt(bc, dpcs * 2);
bigdecimal.DefaultPrecision := dpcs;
// 連立方程式の解法 解法できない場合あり
j := LUecomposition;
if (j = 0) or (j = 2) then exit; // 対角値 0 なら終了
// 解法出来なかったら
if j = 1 then begin
MessageDlg('解法時対角ゼロか、ゼロに近い値が発生しました、入力値を見直して下さい。' + #13#10
+ '有効桁数が不足している場合もあります。' + #13#10
+ '有効桁数が小さい場合は、許容誤差の値を見直して下さい。', mtInformation,[mbOk], 0);
exit;
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_Bigdecimal.zip
三角関数、逆三角関数、その他関数、 連立一次方程式の解法 に戻る