Форум для моих друзей
http://atheist4.flybb.ru/

Нахождение концентрации вещества при диффузии.
http://atheist4.flybb.ru/topic31.html		 Страница 1 из 1
Автор:  atheist4 [ 02-06, 08:26 ]
Заголовок сообщения:  Нахождение концентрации вещества при диффузии.
Обновлены публикации:

http://atheist4.narod.ru/diffuz.htm

http://atheist4.narod.ru/simplediffuz.htm

Исправлена опечатка.

Хотел решение трансцендентного уравнения искать методом касательных, а не методом деления отрезка пополам. Но погрешность вычислений стала больше, и программа зависает уже при N=1000, при h=1, L=3, DX=1E-8 ,например:

procedure PCoInvoluta;
const
DX=1E-8;
var
h:real;
L:real;
M:real;
K:integer;
W: real;
N: integer;

function Y(X:real):real;
begin
Y:=h*Cos(X*L)/Sin(X*L)-X;
end;

function Z(X:real):real;
begin
Z:=-h*L/Sin(X*L)/Sin(X*L)-1;
end;

begin
h:= StrToFloat(Form1.Edit1.Text);
L:= StrToFloat(Form1.Edit2.Text);
N:= StrToInt(Form1.Edit3.Text);
if Form1.CheckBox1.State = cbChecked
then Form1.Memo1.Visible:=true
else Form1.Memo1.Visible:=false;
Form1.Memo2.Visible:=false;
Form1.Memo3.Visible:=false;
for K:=1 to N do
begin
M:=Pi/4/L;
if Y(Pi*(K-1)/L+M)<=0
then
repeat M:=M/2;
until Y(Pi*(K-1)/L+M)>0;
W:=Pi*(K-1)/L+M;
While Y(W)/(h*L+1)>DX do W:=W-Y(W)/Z(W);
Form1.Memo1.Lines.Add('W('+IntToStr(K)+')='+FloatToStrF(W,ffFixed,20,14));
Form1.Memo2.Lines.Add('d('+IntToStr(K)+')='+FloatToStr(Y(W)));
Form1.Memo3.Lines.Add('f('+IntToStr(K)+')='+FloatToStr(abs(Y(W)/W)));
end;
Form1.Memo1.Lines.Add('Êîðíè óðàâíåíèÿ W=h*ctg(W*L)');
Form1.Memo1.Lines.Add('h='+FloatToStr(h)+' L='+FloatToStr(L));
Form1.Memo1.Visible:=true;
end;

Поэтому оставлю лучше всё, как было - буду искать методом деления отрезка пополам:

procedure PCoInvoluta;
const
DX=1E-30;{погрешность вычисления}
var
h:real;{Значение h}
L:real;{Значение L}
M:real;
K:integer;
I,P: integer;
X1,X2 :real;
W: real;
N: integer;
F: TextFile;
function Y(Z:real):real;
begin
Y:=h*Cos(Z*L)/Sin(Z*L)-Z;
end;
begin
h:= StrToFloat(Form1.Edit1.Text);
L:= StrToFloat(Form1.Edit2.Text);
N:= StrToInt(Form1.Edit3.Text);
if Form1.CheckBox1.State = cbChecked
then Form1.Memo1.Visible:=true
else Form1.Memo1.Visible:=false;
Form1.Memo2.Visible:=false;
Form1.Memo3.Visible:=false;
Form1.Memo1.Lines.Add('Корни уравнения W=h*ctg(W*L)');
Form1.Memo1.Lines.Add('h='+FloatToStr(h)+' L='+FloatToStr(L));
Form1.Memo2.Lines.Add('Абсолютные погрешности:');
Form1.Memo3.Lines.Add('Относитенльные погрешности:');
if Form1.CheckBox2.State = cbChecked
then Form1.Label5.Caption:='h'+Form1.Edit1.Text+'L'+Form1.Edit2.Text+'.txt';
for K:=1 to N do
begin
M:=Pi/4/L;
if Y(Pi*(K-1)/L+M)<=0
then
repeat M:=M/2;
until Y(Pi*(K-1)/L+M)>0;
X1:=Pi*(K-1)/L+M;
X2:=Pi*(K-1)/L+2*M;
W:=X2;
P:=Round(Ln(M/DX)/Ln(2));
for I:=1 to P do
begin
if Y(X1)*Y(W)<=0
then X2:=W
else X1:=W;
W:=(X1+X2)/2;
end;
Form1.Memo1.Lines.Add('W('+IntToStr(K)+')='+FloatToStr(W));
Form1.Memo2.Lines.Add('d('+IntToStr(K)+')='+FloatToStr(Y(W)));
Form1.Memo3.Lines.Add('f('+IntToStr(K)+')='+FloatToStr(abs(Y(W)/W)));
end;
Form1.Memo1.Visible:=true;
if Form1.CheckBox2.State = cbChecked
then
begin
AssignFile(F,Form1.Label5.Caption);
Rewrite(F);
WriteLn(F,Form1.Memo1.Text);
WriteLn(F,Form1.Memo2.Text);
WriteLn(F,Form1.Memo3.Text);
Close(F);
end;
end;
Автор:  atheist4 [ 02-06, 22:26 ]
Заголовок сообщения: 
Обновлённая программа

http://atheist4.narod.ru/programs/diffuz.rar

показывает терерь проверку решения дифференциального уравнения

procedure PDiffuz;
const
DX=1E-17; { погрешность вычисления }
var
h: real; { Значение h }
L: real; { Значение L - длина промежутка }
b: real; { Значение b }
D: real; { Значение коэффициента диффузии D }
A0: real; { Значение начальной концентрации }
N: integer; { Количество членов ряда}
M: real;
K: integer; { Номер члена ряда }
I, P: integer;
X1, X2: real;
x: real; { Координата x }
t:real; { Время t }
W: real; { Корни уравнения W=h*ctg(W*L) }
a: real; { Коэффициенты разложения единицы в ряд по cos(wx) }
V: real; { K-й член ряда в сумме n(x, t) }
Q: real; { n(x, t) - концентрация }
G: real; { K-й член ряда в сумме dn(x, t)/dx }
U: real; { dn(x, t)/dx - производная от концентрации по координате x }
{ Данные для проверки правильности решения диф. уравнения }
C1: real; { K-й член ряда в сумме dn(x, t)/dt }
F1: real; { dn(x, t)/dt - производная от концентрации по времени t}
C2: real; { K-й член ряда в сумме d2n(x, t)/dx2 }
F2: real; { d2n(x, t)/dx2 - вторая производная от концентрации по x }

function Y(Z:real):real;
begin
Y:=h*Cos(Z*L)/Sin(Z*L)-Z;
end;

begin
U:=0;
L:= StrToFloat(Form1.Edit1.Text);
h:= StrToFloat(Form1.Edit2.Text);
b:= StrToFloat(Form1.Edit3.Text);
D:= StrToFloat(Form1.Edit4.Text);
A0:=StrToFloat(Form1.Edit5.Text);
x:= StrToFloat(Form1.Edit6.Text);
t:= StrToFloat(Form1.Edit7.Text);
N:= StrToInt(Form1.Edit8.Text);
Q:= A0*exp(-b*t);
U:= 0;
F1:= -A0*b*exp(-b*t);
F2:= 0;
if x>L then
begin
x:=L;
Form1.Edit6.Text:=FloatToStr(x);
end;
for K:=1 to N do
begin
M:=Pi/4/L;
if Y(Pi*(K-1)/L+M)<0>0;
X1:= Pi*(K-1)/L+M;
X2:= Pi*(K-1)/L+2*M;
W:= X2;
P:= Round(Ln(M/DX)/Ln(2));
for I:=1 to P do
begin
if Y(X1)*Y(W)<0>L then
begin
x:=L;
Form1.Edit6.Text:=FloatToStr(x);
end;
for K:=1 to N do
begin
W:= Pi*(2*K-1)/2/L;
a:= sin(W*L)/(W*L/2+sin(2*W*L)/4);
V:= A0*b*a/(D*W*W-b)*(exp(-b*t)-exp(-D*W*W*t))*cos(W*x);
Q:= Q+V;
G:= -A0*b*a*W/(D*W*W-b)*(exp(-b*t)-exp(-D*W*W*t))*sin(W*x);
U:= U+G;
C1:= A0*b*a/(D*W*W-b)*(-b*exp(-b*t)+D*W*W*exp(-D*W*W*t))*cos(W*x);
F1:= F1+C1;
C2:= -A0*b*a*W*W/(D*W*W-b)*(exp(-b*t)-exp(-D*W*W*t))*cos(W*x);
F2:= F2+C2;
end;
Form1.Label10.Caption:=' Концентрация n(x, t) = '+FloatTostr(Q)+
#13+#13+'Градиент концентрации dn(x,t)/dx = '+ FloatTostr(U) ;
Form1.Label11.Caption:='Проверка: dn/dt = '+ FloatToStr(F1)+ #13+
' D*d2n/dx2 = '+ FloatToStr(D*F2);
if x=L then
Form1.Label11.Caption:=Form1.Label11.Caption + #13 +
'на границе A*exp(-b*t) = ' +FloatToStr(A0*exp(-b*t));

end;

Эта программа даёт верное решение дифференциального уравнения.

http://atheist4.narod.ru/diffuz.htm
Автор:  atheist4 [ 02-06, 22:34 ]
Заголовок сообщения: 
Но такая же программа с проверкой для второй задачи даёт ошибочное решение дифференциального уравнения

procedure PSimpleDiffuz;
var
h: real; { Значение h }
L: real; { Значение L - длина промежутка }
b: real; { Значение b }
D: real; { Значение коэффициента диффузии D }
A0: real; { Значение начальной концентрации }
N: integer; { Количество членов ряда}
K: integer; { Номер члена ряда }
x: real; { Координата x }
t:real; { Время t }
W: real; { Корни уравнения cos(W*L)=0 }
a: real; { Коэффициенты разложения единицы в ряд по cos(wx) }
V: real; { K-й член ряда в сумме n(x, t) }
Q: real; { n(x, t) - концентрация }
G: real; { K-й член ряда в сумме dn(x, t)/dx }
U: real; { dn(x, t)/dx - производная от концентрации по координате x }
{ Данные для проверки правильности решения диф. уравнения }
C1: real; { K-й член ряда в сумме dn(x, t)/dt }
F1: real; { dn(x, t)/dt - производная от концентрации по времени t}
C2: real; { K-й член ряда в сумме d2n(x, t)/dx2 }
F2: real; { d2n(x, t)/dx2 - вторая производная от концентрации по x }

function Y(Z:real):real;
begin
Y:=h*Cos(Z*L)/Sin(Z*L)-Z;
end;

begin
U:=0;
L:= StrToFloat(Form1.Edit1.Text);
h:= StrToFloat(Form1.Edit2.Text);
b:= StrToFloat(Form1.Edit3.Text);
D:= StrToFloat(Form1.Edit4.Text);
A0:=StrToFloat(Form1.Edit5.Text);
x:= StrToFloat(Form1.Edit6.Text);
t:= StrToFloat(Form1.Edit7.Text);
N:= StrToInt(Form1.Edit8.Text);
Q:= A0*exp(-b*t);
U:= 0;
F1:= -A0*b*exp(-b*t);
F2:= 0;
if x>L then
begin
x:=L;
Form1.Edit6.Text:=FloatToStr(x);
end;
for K:=1 to N do
begin
W:= Pi*(2*K-1)/2/L;
a:= sin(W*L)/(W*L/2+sin(2*W*L)/4);
V:= A0*b*a/(D*W*W-b)*(exp(-b*t)-exp(-D*W*W*t))*cos(W*x);
Q:= Q+V;
G:= -A0*b*a*W/(D*W*W-b)*(exp(-b*t)-exp(-D*W*W*t))*sin(W*x);
U:= U+G;
C1:= A0*b*a/(D*W*W-b)*(-b*exp(-b*t)+D*W*W*exp(-D*W*W*t))*cos(W*x);
F1:= F1+C1;
C2:= -A0*b*a*W*W/(D*W*W-b)*(exp(-b*t)-exp(-D*W*W*t))*cos(W*x);
F2:= F2+C2;
end;
Form1.Label10.Caption:=' Концентрация n(x, t) = '+FloatTostr(Q)+
#13+#13+'Градиент концентрации dn(x,t)/dx = '+ FloatTostr(U) ;
Form1.Label11.Caption:='Проверка: dn/dt = '+ FloatToStr(F1)+ #13+
' D*d2n/dx2 = '+ FloatToStr(D*F2);
if x=L then
Form1.Label11.Caption:=Form1.Label11.Caption + #13 +
'на границе A*exp(-b*t) = ' +FloatToStr(A0*exp(-b*t));

end;

Где ошибка в решении второй задачи:

http://atheist4.narod.ru/simplediffuz.htm ?
Автор:  atheist4 [ 09-06, 20:00 ]
Заголовок сообщения: 
Всё проверено. Ошибки нет. Просто на границе ряд не допускает почленного дифференцирования. Обновлённые программы можете скачать по ссылкам:
http://atheist4.narod.ru/programs/diffuz.rar

http://atheist4.narod.ru/programs/simplediffuz.rar

http://atheist4.narod.ru/programs/diffuzall.rar
Страница 1 из 1 Часовой пояс: UTC + 3 часа
Powered by phpBB® Forum Software © phpBB Group
http://www.phpbb.com/