Расчёт коэффицента диффузии D
есть уравнения вида:
C1(x)=C0*erfc(x1/2sqrt(D*t1))
C2(x)=C0*erfc(x2/2sqrt(D*t2))
Алгоритм вычисления:
1.Задается время отжига и начальное значение С0 в первом цикле(к=1)
2.Ввод зависимости С(х) из n точек
3.Вычисления отношения Сn(x)/C0
4.Вычисление аргумента Z функции erfc(Z)=Cn/C0 (xn/2sqrt(D*t) = Z)
5.Dn = xn*xn/(4*t*Z*Z)
6.Вычисление среднего Dk = сумме по n (Dn / n)
7. |Dk - D(k-1)|<Dk/100
нет да
k:=k+1 конец программы
Co(k+1):= C1(x)/erfc(x1/2sqrt(Dk * t))
Вот код:
unit Unit1;
interface
uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, StdCtrls, Math,
ap, normaldistr;
type
TForm1 = class(TForm)
Edit1: TEdit;
Edit2: TEdit;
Edit4: TEdit;
Edit6: TEdit;
Label1: TLabel;
Label2: TLabel;
Label3: TLabel;
Label6: TLabel;
Button1: TButton;
Edit7: TEdit;
Label8: TLabel;
Label9: TLabel;
Label10: TLabel;
Label11: TLabel;
Label12: TLabel;
Label13: TLabel;
Label14: TLabel;
Edit5: TEdit;
Edit3: TEdit;
Label4: TLabel;
procedure Edit7MouseUp(Sender: TObject; Button: TMouseButton;
Shift: TShiftState; X, Y: Integer);
procedure Button1Click(Sender: TObject);
procedure FormCreate(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;
const
n=20;
var
Form1: TForm1;
k : integer;
// no : integer;
Co,t,Dks : real;
dd,D,d1,d2 : real;
C1,C2 : real;
x1,x2 : real;
z1,z2 : real;
// i: integer;
{Co,t,Dks : Extended;
dd,D,d1,d2 : Extended;
C1,C2 : Extended;
x1,x2 : Extended;
z1,z2 : Extended;}
implementation
{$R *.dfm}
procedure TForm1.Button1Click(Sender: TObject);
begin
//Начальные данные РАссматриваем вариант 2 точек!
Co := StrToFloat(Edit1.Text);
C1 := StrToFloat(Edit2.Text);
x1 := StrToFloat(Edit4.Text);
C2 := StrToFloat(Edit3.Text);
x2 := StrToFloat(Edit5.Text);
t := StrToFloat(Edit6.Text);
k:=0;
// no := StrToInt(Edit7.Text);
Dks:=0;
repeat
if Co=0 then
begin
Label4.visible := true; break;
Label4.Caption := 'Co обратилось в ноль! Плохо!';
end
else
begin
//для каждого i считается считается Zi
z1:= 1/2*sqrt(Pi)*( C1/Co + Pi*Power(C1/Co,3)/12 + 7*Pi*Pi*Power(C1/Co,5)/480 +
127*Power(Pi,3)*Power(C1/Co,7)/40320 + 4369*Power(Pi,4)*Power(C1/Co,9)/5806080 +
34807*Power(Pi,5)*Power(C1/Co,11)/182476800 );
z2:= 1/2*sqrt(Pi)*( C2/Co + Pi*Power(C2/Co,3)/12 + 7*Pi*Pi*Power(C2/Co,5)/480 +
127*Power(Pi,3)*Power(C2/Co,7)/40320 + 4369*Power(Pi,4)*Power(C2/Co,9)/5806080 +
34807*Power(Pi,5)*Power(C2/Co,11)/182476800 );
end;
// сумма по i . у меня оно равно 2.
d1 := (x1*x1) / (4*t*z1*z1);
d2 := (x2*x2) / (4*t*z2*z2);
D:= (d1+d2)/2; dd := sqrt(d1*d1+d2*d2)/2;
if abs(D-Dks) < D/100 then
begin//конец программы:)
break;
end
else
begin
k:=k+1;
Dks:=D;
Co:= C1/erfc(x1 / (2 * sqrt(t*d)));
end;
until k=20000;
Label4.Caption := FloatToStr(Co); label12.Caption := FloatToStr(D);
label13.Caption := FloatToStr(dd); label14.Caption := IntToStr(K);
label8.Visible:=true; label9.Visible:=true; label10.Visible:=true; label11.Visible:=true;
label12.Visible:=true; label13.Visible:=true; label14.Visible:=true; label4.Visible:=true;
end;
procedure TForm1.Edit7MouseUp(Sender: TObject; Button: TMouseButton;
Shift: TShiftState; X, Y: Integer);
begin
Edit7.Text := ' ';
end;
procedure TForm1.FormCreate(Sender: TObject);
begin
Edit7.visible := false;
Edit1.Text := FloatToStr(33); //Co
Edit2.Text := FloatToStr(16.3); //C1
Edit3.Text := FloatToStr(8.5); //C2
Edit4.Text := FloatToStr(0.4); //X1
Edit5.Text := FloatToStr(0.8); //X2
Edit6.Text := IntToStr(334800); //t
end;
end.
на всякий случай что такое erfc :
http://ru.wikipedia.org/wiki/Функция_ошибок
