velocity computing.pas

二分法解非线性方程－－－相变过程中的界面迁移速率

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, ExtCtrls, ComCtrls;

type
  TFloatArray = array of extended;
  TForm1 = class(TForm)
    Button1: TButton;
    LabeledEdit1: TLabeledEdit;
    LabeledEdit2: TLabeledEdit;
    LabeledEdit3: TLabeledEdit;
    LabeledEdit4: TLabeledEdit;
    LabeledEdit5: TLabeledEdit;
    LabeledEdit6: TLabeledEdit;
    LabeledEdit7: TLabeledEdit;
    LabeledEdit8: TLabeledEdit;
    LabeledEdit9: TLabeledEdit;
    LabeledEdit10: TLabeledEdit;
    LabeledEdit11: TLabeledEdit;
    LabeledEdit12: TLabeledEdit;
    LabeledEdit13: TLabeledEdit;
    LabeledEdit14: TLabeledEdit;
    Memo1: TMemo;
    CheckBox1: TCheckBox;
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
    procedure DDHRT(A,B,H,EPS:extended;var X:TFloatArray;N:integer;var M:integer);
    function  F(V:extended):extended;
  end;

var
  Form1: TForm1;
  DeltaT,DeltaH,R,Kl,K,m0,D,C00,Gama,r1,Tm:extended;

implementation

{$R *.dfm}

procedure TForm1.DDHRT(A,B,H,EPS:extended;var X:TFloatArray;N:integer;var M:integer);
label
  T10,T20;
var
  Z,Y,Z1,Y1,Z0,Y0:extended;
begin
 M:=0;
 Z:=A;
 Y:=F(Z);
T10: if((Z>(B+H/2.0))or(M=N))then
   exit;
 if(abs(Y)<EPS)then
 begin
   inc(M);
   X[M-1]:=Z;
   Z:=Z+H/2.0;
   Y:=F(Z);
   goto T10;
 end;
 Z1:=Z+H;
 Y1:=F(Z1);
 if(abs(Y1)<EPS)then
 begin
   inc(M);
   X[M-1]:=Z1;
   Z:=Z1+H/2.0;
   Y:=F(Z);
   goto T10;
 end;
 if(Y*Y1>0.0)then
 begin
   Y:=Y1;
   Z:=Z1;
   goto T10;
 end;
 T20: if(abs(Z1-Z)<EPS)then
 begin
   inc(M);
   X[M-1]:=(Z1+Z)/2.0;
   Z:=Z1+H/2.0;
   Y:=F(Z);
   goto T10;
 end;
 Z0:=(Z1+Z)/2.0;
 Y0:=F(Z0);
 if(abs(Y0)<EPS)then
 begin
   inc(M);
   X[M-1]:=Z0;
   Z:=Z0+H/2.0;
   Y:=F(Z);
   goto T10;
 end;
 if(Y*Y0<0.0)then
 begin
   Z1:=Z0;
   Y1:=Y0;
 end
 else
 begin
   Z:=Z0;
   Y:=Y0;
 end;
 goto T20;
end;

function TForm1.F(V:extended):extended;
begin
  Result:=V*DeltaH*R/Kl-DeltaT+m0*(K-1)*R*V/D*C00/(1-(1-K)*(4/(3+sqrt(1+2*D*r1/sqr(V)))))+2*Gama*Tm/R;
//  Result:=(((((V-5)*V+3)*V+1)*V-7)*V+7)*V-20.0;
//  Result:=exp(V)-3.0*V;
end;

procedure TForm1.Button1Click(Sender: TObject);
const
  N=20;
  EPS=1.0e-12;
var
  X:TFloatArray;
  A,B,H:extended;
  M,I:integer;
begin
 SetLength(X,N);
 for I:=0 to High(X)do
   X[I]:=0.0;
// DDHRT(-2,5,0.2,1.0e-6,X,6,M);
// DDHRT(0,1,0.0001,1.0e-5,X,50,M);
 DeltaT:=StrToFloat(LabeledEdit1.Text);
 DeltaH:=StrToFloat(LabeledEdit2.Text);
 R:=StrToFloat(LabeledEdit3.Text);
 Kl:=StrToFloat(LabeledEdit4.Text);
 K:=StrToFloat(LabeledEdit5.Text);
 m0:=StrToFloat(LabeledEdit6.Text);
 D:=StrToFloat(LabeledEdit7.Text);
 C00:=StrToFloat(LabeledEdit8.Text);
 Gama:=StrToFloat(LabeledEdit9.Text);
 r1:=StrToFloat(LabeledEdit10.Text);
 Tm:=StrToFloat(LabeledEdit11.Text);

 M:=0;
 SetLength(X,N);
 for I:=0 to High(X)do
   X[I]:=0.0;
 A:=StrToFloat(LabeledEdit12.Text);
 B:=StrToFloat(LabeledEdit13.Text);
 H:=StrToFloat(LabeledEdit14.Text);
 DDHRT(A,B,H,EPS,X,N,M);//}
 for I:=0 to M-1 do
 begin
//   ListBox1.Items.Add(FloatToStr(I+1)+'  '+FloatToStr(X[I]));
   Memo1.Lines.Add(FloatToStr(I+1)+#9+FloatToStr(X[I]))
 end;
 if CheckBox1.Checked then
   showmessage('over!');
end;

end.