連立一次方程式の解法8

 連立一次方程式の解法の続きです。
今回は、Gauss Jordan法 です。
 この計算は、各種計算プログラム例の中のガウスの消去法と同じものです、条件シンボル {$DEFINE NNM2} {$DEFINE NNM1} と変更すれば、Gause Jordan法と同じ計算になります。
但し、対角の値(ピボットの値)の判定値が固定となっている為、判定ミスの可能性があります。
 此処でのプログラムは、検算による判定として、判定ミスの防止をしています。
epsilonの値と、行列aの値から対角値(ピボット)の判定値を設定しても良いのですが、設定が難しいので採用していません。

次の例は4×6の行列入力例です。
  8, 16, 2432, 160, 1
  27, 1217,   70, 2
  6, 17, 3259, 198, 3
  7, 22, 46, 105, 291, 4
黒字は、本来の4×4の行列でaの値、赤文字はbの値、青文字は行の番号です。

演算結果
 100041
 010034
 0010,  22
 000113
黒文字の行列は、対角の値が全て1になり、赤文字の列には、Xの値が入ります。
青文字の部分は、行の番号ですが、これでピボット操作により行がどの様に移動したか確認が出来ますが、演算とは直接関係が無いので省略可能です。

プログラム

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;
    LabeledEdit2: TLabeledEdit;
    CheckBox1: TCheckBox;
    procedure BitBtn1Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure BitBtn2Click(Sender: TObject);
    procedure RadioGroup1Click(Sender: TObject);
  private
    { Private 宣言 }
    function Gauss_Jordan: 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 = ((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 : array of array of extended;  // A行列配列
  a    : array of array of extended;
  b    : array of extended;           // b行列

  allowable_errorr : extended;
  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;

// ガウス・ジョルダンの消去法による連立方程式の解法
function TForm1.Gauss_Jordan: byte;
const
  sp = ' ';
var
  i, j, k, n, pivot, h : integer;
  sum, amax, p, d : extended;
  x, y: array of extended;
  LL, spl : string;
begin
  setlength(a, Nl, Nl + 2);                   // 前進消去計算用
  setlength(x, Nl);                           // 答えx値
  setlength(y, Nl);                           // 検算差分値

  memo1.Clear;

  // matA,bデーターを計算用にコピー
  for j := 0 to Nl - 1 do begin
    for k := 0 to Nl - 1 do
      a[j, k] := matA[j, k];
    a[j, Nl] := b[j];                         // aのNl列のにbの値セット
    a[j, Nl + 1] := j + 1;                    // Nl+1列matA行位置初期化
  end;

  p := 0;                                        // 単なる初期化

  // 前進消去 Gauss Jordan法
  for i := 0 to Nl - 1 do begin
    pivot := i;
    // pivot行から最大値検索                     // pivot行
    amax := abs(a[pivot, pivot]);
    for k := pivot + 1 to Nl - 1 do begin
      if abs(a[k, i]) > amax then begin
        pivot := k;
        amax := abs(a[pivot, i]);
      end;
    end;
    // 最大値行が違ったら行交換
    if pivot <> i then begin
      for k := 0 to Nl + 1 do begin
        sum := a[i, k];
        a[i, k] := a[pivot, k];
        a[pivot, k] := sum;
      end;
    end;

    // 対角の値がゼロだったらブレーク
    p := a[i, i];                                // 対角の値
    if p = 0 then break;

    // 対角の値でその行の除算
    for j := 0 to Nl do a[i, j] := a[i, j] / p;  // bの値も含む

    // 対角化
    for j := 0 to Nl - 1 do begin
      if i <> j then begin
        d := a[j, i];
        for k := i to Nl do a[j, k] := a[j, k] - d * a[i, k];
      end;
    end;
  end;

  // 前進消去中対角にゼロがでたら演算中止
  if p = 0 then begin
    result := 1;                                 // 戻り値 1
    exit;
  end;

  // xの値取り出し  Nl列がxの値
  for i := 0 to Nl - 1 do x[i] := a[i, Nl];

  // 各行X値代入答えの差分計算
  for j := 0 to Nl - 1 do begin
    sum := 0;
    for i := 0 to Nl - 1 do
      sum := sum + matA[j, i] * x[i];          // 行の値
    y[j] := b[j] - sum;                        // 差分
  end;

  // 演算中の行の入れ替えによる行の位置表示  Nl+1の列が行位置
  LL := '計算順 ';
  for i := 0 to Nl - 1 do begin
    if i < NL - 1 then LL := LL + floatTostr(a[i, Nl + 1]) + ','
                  else LL := LL + floatTostr(a[i, Nl + 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[k]));
  end;

  // 差分判定
  h := 0;
  for k := 0 to Nl - 1 do begin
    if abs(b[k]) > allowable_errorr then sum := y[k] / b[k]
                                    else sum := y[k];
    if abs(sum) > allowable_errorr then inc(h);
  end;

  if checkbox1.Checked = false then begin
    if h = 0 then memo1.Lines.Append('判定   OK')
             else memo1.Lines.Append('判定   NG');
    if h > 0 then begin
      result := 3;                               // 戻り値 3 判定NG 終了
      exit;
    end;
  end;

  // xの値表示
  for k := 0 to Nl - 1 do
    memo1.Lines.Append('X' + inttostr(k + 1) + '=' + floatTostr(x[k]));

  result := 0;                                   // 戻り値 0 正常終了
end;

procedure TForm1.BitBtn1Click(Sender: TObject);
var
  i, j, ch : integer;
  maxb, minb, bc : extended;
  acmax : array of extended;
  acmin : array of extended;
begin
  setlength(matA, Nl, 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], 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;
  // 検算結果の判定値設定
  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);      // 誤差判定値

  // 連立方程式の解法  解法できない場合あり
  j := Gauss_Jordan;

  // 解法出来なかったら
  if j = 1 then
    MessageDlg('解法時対角ゼロが発生しました、入力値を見直して下さい。', mtInformation,[mbOk], 0);
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);
var
  a : extended;
begin
  testset;
  epsilon := 1;
  repeat
    epsilon := epsilon / 2;
    a := 1 - epsilon;
  until a = 1;
  epsilon := epsilon * 2;
end;

end.

download Gauss_Jordan.zip

  三角関数、逆三角関数、その他関数、 連立一次方程式の解法 に戻る