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

この発言に対し以下のコメントが寄せられています

#00351　Fermion         さん　RE^2:数値積分用Unit






　積分計算用のUnitを作りました．台形近似，ルジャンドル・ガウス積分，シン
プソンの公式に対応しています．ルジャンドル・ガウス積分，シンプソンの公式
は，Fermionさんのサンプルを引用いたしました．Fermionさん，ありがとうござ
います．

●説明
　ある被積分関数f(x)の，区間[a,b]における定積分の近似値を計算します．台
形近似，シンプソンの公式では，区間[a,b]を何分割して計算するか設定できま
す．ルジャンドル・ガウス積分は，4点公式です．

●注意
　例外発生に対する処理は省略してあります．

●使い方

[例]

implementation

uses UnitIntegral;  //積分計算用のUnitを参照できるようにします

{$R *.DFM}

function SomeFuncion(x: Double) : Double; //任意の関数を設定する
begin
  Result := 1/(1+x*x);//例えば
end;


//=====メインルーチン======台形近似の場合
procedure TForm1.Button1Click(Sender: TObject);
var
  Integral: TIntegral; // TIntegral型のインスタンス
  a,b: Double;
  n: Integer;
begin
  //初期設定
  a := 1;//例えば
  b := 2;//例えば
  n := 100;//例えば

  Integral := TIntegral.Create(a,b,n,SomeFunction);//変数をセット
  Edit1.Text := FloatToStr(Integral.Show);//Showメソッドで，
                                          //計算結果を呼び出す
  Integral.Free;//解放
end;


//=====メインルーチン======ルジャンドル・ガウスの場合
procedure TForm1.Button1Click(Sender: TObject);
var
  Integral: TLGIntegral; // TLGIntegral型のインスタンス
  a,b: Double;
begin
  //初期設定
  a := 1;//例えば
  b := 2;//例えば

  Integral := TLGIntegral.Create(a,b,SomeFunction);//変数をセット
  Edit1.Text := FloatToStr(Integral.Show);//Showメソッドで，
                                          //計算結果を呼び出す
  Integral.Free;//解放
end;


//=====メインルーチン======シンプソンの公式の場合
procedure TForm1.Button1Click(Sender: TObject);
var
  Integral: TSimpIntegral; // TSimpIntegral型のインスタンス
  a,b: Double;
  n: Integer;
begin
  //初期設定
  a := 1;//例えば
  b := 2;//例えば
  n := 100;//例えば

  Integral := TSimpIntegral.Create(a,b,n,SomeFunction);//変数をセット
  Edit1.Text := FloatToStr(Integral.Show);//Showメソッドで，
                                          //計算結果を呼び出す
  Integral.Free;//解放
end;



//##########################################積分計算Unitのコード
unit UnitIntegral;

interface

type
  TFunction = function(x: Double) : Double; // 関数型の宣言
//注意：コンパイル時，メインルーチン中のCreateメソッドで，
//「変数に互換性がない」とエラーが出たら，
//この部分を以下のようにすると良いでしょう．
//TFunction = function(x: Double) : Double of Object; // 関数型の宣言
  TIntegral = class // class型の宣言 台形近似
  private // privete宣言
    Fa: Double; // 定義域　下
    Fb: Double; // 定義域　上
    Fn: Integer; // 定義域をn等分
    Ff: TFunction; // 関数型のインスタンス
    function Calculate: Double; virtual;// 計算部分
  public // public宣言
    constructor Create(a,b: Double; n: Integer; f: TFunction);
    function Show: Double; //結果の表示
  end;

  TLGIntegral = class(TIntegral) //ルジャンドル・ガウス
  private
    function Calculate: Double; override;//計算部分
  public
    constructor Create(a,b: Double; f: TFunction); //4点
  end;

  TSimpIntegral = class(TIntegral)//シンプソン
  private
    function Calculate: Double; override;//計算部分
  end;

implementation

//####################################################台形近似
constructor TIntegral.Create(a,b: Double; n: Integer; f: TFunction);
begin
  Fa := a; // fieldに代入
  Fb := b;
  Ff := f;
  Fn := n;
end;


function TIntegral.Calculate : Double; // 台形近似計算
var
  i: Integer;
  S,
  xi,xi_1,xdelta: Double;
begin
  S:=0;

  // xの増分幅xdeltaを求める
  xdelta := ( Fb - Fa )/Fn;

  For i:= 0 to Fn - 1 Do // 面積計算
  begin
    xi := Fa + xdelta*i;
    xi_1 := xi + xdelta;
    S := S +( Ff(xi) + Ff(xi_1) )*xdelta/2
  end;

  Result := S;
end;

function TIntegral.Show: Double; // 計算結果を表示
begin
  Result := Calculate;
end;


//####################################################ルジャンドル・ガウ
ス
constructor TLGIntegral.Create(a,b: Double; f: TFunction);
begin
  Fa := a;
  Fb := b;
  Ff := f;
end;

function TLGIntegral.Calculate: Double;
const
  T : array[ 1..4 ] of Double =   ( -0.861136,  // T1
                                    -0.339981,  // T2
                                     0.339981,  // T3
                                     0.861136 );// T4
  C : array[ 1..4 ] of Double =   (  0.347855,  // C1
                                     0.652145,  // C2
                                     0.652145,  // C3
                                     0.347855 );// C4
var
  BmA, ApB, SIGMA, X : Double;
  i : Integer;
begin
  BmA := ( Fb - Fa ) / 2.0;
  ApB := ( Fa + Fb ) / 2.0;
  SIGMA := 0;
  for i := 1 to 4 do
    begin
      X := ApB + BmA * T[i];
      SIGMA := SIGMA + C[i] * Ff( X );
    end;
  Result := BmA * SIGMA;
end;


//####################################################シンプソン
function TSimpIntegral.Calculate: Double;
var
  h, S1, S2, X : Double;
  i : Integer;
begin
  h := ( Fb - Fa ) / ( 2 * Fn );
  S1 := 0;
  S2 := 0;
  X := Fa;
  for i := 1 to Fn - 1 do
    begin
      X := X + h;
      S1 := S1 + Ff( X );
      X := X + h;
      S2 := S2 + Ff( X );
    end;
  S1 := S1 + Ff( X + h );
  Result := ( Ff( Fa ) + 4.0 * S1 + 2.0 * S2 + Ff( Fb ) )
              * h / 3.0;
end;

end.


--------------------------------------------------------------
_/_/_/              98/01/16(金)  pm 09:23              _/_/_/
_/_/_/      はんばあぐ(BZQ00131@niftyserve.or.jp)       _/_/_/

Original document by はんばあぐ      氏 ID:(BZQ00131)

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