この発言は
#00218　Fermion         さんの数値積分 << シンプソン >>
に対するコメントです






■説明
　ルジャンドル・ガウス法、シンプソン法、台形近似法による数値積分
を数値積分ユニットにまとめました。サンプルコードの“関数ユニット”
部分を別ユニットにして下さい。

  [a, b] に関する積分値Sを求める場合、

  S := Integral( 関数, 方式, a, b, 分割数 );

　とします。なお、方式の指定は列挙型

  type
    TFormula   = ( fLeGauss, fSimpson, fTrapezoidal );

　を用いて下さい。{S:=Integral( 関数, fSimpson, a, b, 分割数 );}


■使用例
  Form1 に Edit1 〜 Edit4, RadioGroup1, Button1 を適当に配置して
下さい。
　　　 b
　Ｓ＝∫f(X)dX  について、
       a
　　Edit1 : a 入力用、  Edit2 : b 入力用
　　Edit3 : n 入力用（分割数。Simpson's rule, Trapezoidal rule ）
　　Edit4 : Ｓ（積分結果）表示用

    RadioGroup1.Items[0] → 'Legendre Gauss'
    RadioGroup1.Items[1] → 'Simpson'
    RadioGroup1.Items[2] → 'Trapezoidal'

    Button1 : 計算処理開始用

　と設定して、Button1 の OnClickイベントを下記サンプルコードの様
にして下さい。


■注意点
　例外発生、計算結果の有効桁等については検討の必要があります。
　ルジャンドル・ガウスは４点公式を使用しています。


■サンプルコード
//== 使用例 ===========================================================
implementation

uses MathEx; //←※注意

{$R *.DFM}

//*** 非積分関数の定義 ***
function F( X : Extended ): Extended;
begin
  Result := 1 /( 1 + X*X );
end;
//*** 数値積分実行 *******
procedure TForm1.Button1Click(Sender: TObject);
var
  a, b : Extended; //積分区間
  n    : Integer;  //分割数
begin
  a := StrToFloat( Edit1.Text );
  b := StrToFloat( Edit2.Text ); n := StrToInt( Edit3.Text );
  Edit4.Text
   := FloatToStr(Integral(F,TFormula(RadioGroup1.ItemIndex),a,b,n)); 
end;
//== 関数ユニット =====================================================
unit MathEx;

interface

type
  TIntegrand = function( X : Extended ): Extended;
  TFormula   = ( fLeGauss, fSimpson, fTrapezoidal );

  function Integral( Integrand : TIntegrand; Formula : TFormula;
                  const a, b : Extended; const n : Integer ): Extended;

implementation

function Integral( Integrand : TIntegrand; Formula : TFormula;
                  const a, b : Extended; const n : Integer ): Extended;
const
  T : array[ 1..4 ] of Extended = ( -0.861136,  // t1
                                    -0.339981,  // t2
                                     0.339981,  // t3
                                     0.861136 );// t4
  C : array[ 1..4 ] of Extended = (  0.347855,  // c1
                                     0.652145,  // c2
                                     0.652145,  // c3
                                     0.347855 );// c4
var
  BmA, ApB,
  SIGMA1,
  SIGMA2,
  X, h    : Extended;
  i       : Integer;
begin
  Result := 0;
  if not Assigned( Integrand ) then Exit;
  case Formula of
    fLeGauss :     //Legendre Gauss
      begin
        BmA := ( b - a ) / 2.0;
        ApB := ( a + b ) / 2.0;
        SIGMA1 := 0;
        for i := 1 to 4 do
          begin
            X := ApB + BmA * T[ i ];
            SIGMA1 := SIGMA1 + C[ i ] * Integrand( X );
          end;
        Result := BmA * SIGMA1;
      end;
    fSimpson :     //Simpson
      begin
        h := ( b - a ) / ( 2 * n );
        SIGMA1 := 0;
        SIGMA2 := 0;
        X := a;
        for i := 1 to n - 1 do
          begin
            X := X + h;
            SIGMA1 := SIGMA1 + Integrand( X );
            X := X + h;
            SIGMA2 := SIGMA2 + Integrand( X );
          end;
        SIGMA1 := SIGMA1 + Integrand( X + h );
        Result := ( Integrand( a ) + 4.0 * SIGMA1 + 2.0 * SIGMA2
                                     + Integrand( b ) ) * h / 3.0;
      end;
    fTrapezoidal : //Trapezoidal
      begin
        SIGMA1 := 0;
        h := ( b - a ) / n;
        for i := 1 to n - 1 do
          begin
            X := a + i * h;
            SIGMA1 := SIGMA1  + Integrand( X );
          end;
        SIGMA1 := SIGMA1 + ( Integrand( a ) + Integrand( b ) ) / 2.0;
        Result := h * SIGMA1;
      end
    else Result := 0;
  end;
end;
end.
//=====================================================================


                                     98/01/17(土) 17:19 Fermion [KHF03264]

Original document by Fermion         氏 ID:(KHF03264)

ここにあるドキュメントは NIFTY SERVEの Delphi Users' Forum の16番会議室「玉石混淆みんなで作るSample蔵」に投稿されたサンプルです。これらのサンプルはボーランド株式会社がサポートする公式のものではありません。また、必ずしも動作が検証されているものではありません。これらのサンプルを使用したことに起因するいかなる損害も投稿者、およびフォーラムスタッフはその責めを負いません。使用者のリスクの範疇でご使用下さい。