Вопрос

← →
Озадаченный (2002-06-11 15:54) [0]

Уважаемые господа эксперты и просто умные люди!
У меня возник вопрос по типу Extended. Я пишу программу, которая должна найти корень уравнения методом половинного деления в пределах от 1e-30 до 1e30. Если поставить точность меньше чем 1e10 (например 0.1), то программа "виснет", а если точность ставить 1е10, то все работаеть нормально, но точно +-7 миллиардов меня не сильно устраивает. Программа:
...
var a,b,x,n,Fa,Fx:Extended;

...
a:=1e-30;
b:=1e29;
Fa:=формула с параметром а;
repeat
x:=(a+b)/2;
Fx:=формула с параметром х;
if Fa*Fx>0 then a:=x else b:=x;
until (ABS(Fx)<1e5) or (b-a<1e10);
n:=x;
...

Если b-a меньше ставить, то программа виснет. Помогите пожалуйста разобраться.
Спасибо!
С уважением,
Озадаченный

← →
короед (2002-06-11 20:15) [1]

ты уверен, что
if Fa*Fx>0 then a:=x else b:=x;

может
if (Fa-Fx)>0 then a:=x else b:=x;
еще
until (ABS(Fx)<1e5) or ( abs(b-a) <1e10);
вот так. может заработает...

← →
MegaVolt (2002-06-12 14:26) [2]

ИМХО if Fa*Fx>0 then a:=x else b:=x; правильно т.к. эта строка проверяет находится ли точка F(a) и точка F{x} по одну сторону от оси х.
А если прога виснет то есть вероятность что на данном участке есть несколько решений и задавать участок на которосм прога ищет решенние нужно поменьше и с условием что F(а) и F(b) лежат по разные стороны от оси х т.е. одно Ю0 а другое Б0 илши наоборот.
Вот текстик проги которая это всё решает:
function mdop; var x1,x2,f1,f2:real; begin repeat x1:=(a+b-e)/2; x2:=(a+b+e)/2; f1:=f(x1); f2:=f(x2); if f1>f2 then a:=x1 else b:=x2 until abs(a-b)<e*2; mdop:=(a+b)/2; end; end.

← →
Озадаченный (2002-06-13 11:57) [3]

короед © : Да, насчет этого точно уверен, т.к. действительно функция проверяет одинаковые ли знаки получились, а модуль в условии тоже не играет особой роли, т.к. b всегда больше a
MegaVolt © : Я эту функцию и корни вроде бы вычислились, но по-моему они немного неправильные, т.к. я попробовал эту функцию на примере y=x*x-9 и у меня получился корень (в пределах от 0 до 1000) равный 7.9*10 в минус пятой
e:=0.0001; a:=0; b:=1000; repeat x1:=(a+b-e)/2; x2:=(a+b+e)/2; f1:=x1*x1-9; f2:=x2*x2-9; if f1>f2 then a:=x1 else b:=x2; until abs(a-b)<e*2; mdop:=(a+b)/2; showmessage(floattostr(mdop));

а если этим методом, то все нормально:
a:=0; b:=1000; Fa:=a*a-9; repeat x:=(a+b)/2; Fx:=x*x-9; if Fa*Fx>0 then a:=x else b:=x; until (ABS(Fx)<0.1) or (b-a<0.1); showmessage(floattostr(x));

Могу дать формулу: Fx:=x+Na/(2*exp((Ea-Ec-k*T*ln(x/Nc))/(k*T))+1)-ni*ni/x-Nd/(2*exp((Ec+k*T*ln(x/Nc)-Ed)/(k*T))+1);
все переменные известные, кроме переменной x

← →
ZavDim (2002-06-13 12:52) [4]

1. А может решений вообще нет? Тогда ясно, что происходит.
2. Во- вторых мы можем промахнуться с таким диапазоном. Если у Вас будут числа порядка 1е20, то точности представления Extended просто не хватит, чтобы обеспечить Ваши 1е-много.
Вот и смотри, если 1е30, то минимальный доступный разряд 1е10, другие просто не будут присутствовать. Два числа такой величины просто неразличимы в разрядах меньше 10 миллиардов. Вот прога и крутит числа, а нужной разницы не добивается. Если корни порядка 1, то требуемая точность никогда не достигнется. Определись с порядком корней и требуй разумного. Самое лучшее - это не абсолютная, а относительная погрешность.
А про ABS(b-a) - не забывай. Хотя ты от большего отнимаешь меньшее, но все-таки.

← →
MegaVolt (2002-06-14 16:35) [5]

Попробуй пройти по шагам на примере твоей простой функции х*х-9 и посмотри где он неправильно работает но по идее всё должно быть правильно я сам проверял в молодости когда ещё в институте учился так что попробуй и напиши где он начинает неправильно считать

← →
Озадаченный (2002-06-14 18:54) [6]

Вот именно, что он все правильно считает если брать маленькую функцию, даже если пределы брать большие, а в моей функции начинает глючить.

ZavDim: корни по идее должны быть, ведь он же их находит с точностью до 1е10, а дальше не хочет искать... не понимаю. Я делал счет почти по шагам и он "валится" где-то при значении Т=300, видимо, тут для него наступает что-то критическое...
Спасибо большое Вам за ответы, постараюсь учесть все советы и попробовать еще раз.

← →
Озадаченный (2002-06-14 18:58) [7]

я еще не понимаю почему тип Extended не хочет считать до большой точности, хотя может...

← →
ZavDim (2002-06-15 07:13) [8]

Да есть у него (Extended) точность, но относительная!
Попробуй прибавлять к 1 число 1e-39 и посмотри, что будет получаться - так 1 и останется. Он же не резиновый, он хранит ограниченное количество знаков. Если у тебя величины 1e10, то требуемая точность не может быть 1e-20 - хоть лопни (прога просто зациклится).
Возьми любое Extended(sqrt(2)) и посмотри сколько тебе знаков дастся, остальные, ты в этом представлении - под пыткой не вытащищь. Я когда считал тысячу знаков Pi, так использовал другое представление - свое собственное, а иначе ничего не посчитаешь.
x^2-9 - решается, потому, что диапазон валится вниз и приходит к обыденным числам, там можно обеспечить приличную абсолютную точность - но 1e-39 - не верю. 1е-20 еще туда сюда.
Далее. <Не хочет искать...валится> - что значит не хочет: вылетает, ругается - что происходит?

← →
ZavDim (2002-06-15 07:24) [9]

Про точность 1e-39 это я перекрутил, но выложи значения констант в своих функциях - я код покручу. Может и другое что. Ноя слегка по другому ищу.

type
F1=FUNCTION(T:EXTENDED):EXTENDED; {Одноместная функция}

PROCEDURE KORFUN(A,B,EPS:EXTENDED;F:F1;VAR X:EXTENDED;VAR C:BOOLEAN);
VAR B,YA,YB,YX:EXTENDED;
BEGIN
C:=TRUE;
YA:=F(A);
YB:=F(B);
X:=(B+A)/2;
YX:=F(X);
WHILE ABS(YX)>EPS DO
BEGIN
IF ABS(B-A)<SQR(EPSILON) THEN
BEGIN
C:=FALSE;
EXIT;
END;
IF YA*YX<0 THEN
BEGIN
B:=X;
YB:=YX;
END
ELSE
BEGIN
A:=X;
YA:=YX;
END;
X:=(B+A)/2;
YX:=F(X);
END;
END;

При этом, конечно подразумевается, что корень на участке [А,В]есть. Ну, а вызов уж организуешь, я думаю.
Это голый Pascal 7.0.

← →
Озадаченный (2002-06-18 12:54) [10]

Спасибо большое, сейчас попробую, если не получится, то выложу кусок программы с константами, просто я не знаю даже приблизительно в каком диапазоне могут лежать корни, но знаю точно, что больше нуля (ноль нельзя брать, т.к. там деление на ноль происходит), до, по-моему, 10 в 21 степени.
спасибо еще раз за помощь всем откликнувшимся.

← →
Extended (2002-06-18 13:30) [11]

Кстати, в хелпе написано, что Extended 3.6 x 10^–4951 .. 1.1 x 10^4932 так почему же он не вмещает значения хотя бы от 1 до 10^30 с точностью до 1?

← →
Shaman_Naydak (2002-06-19 09:43) [12]

> Extended
Потому, что это не целочисленный тип, а с плавающей запятой..
И точность представления - это 19-20 знаков, а уж сколько нуликов что в одну, что в другую сторону - це дело хозяйское

← →
Lord Warlock (2002-06-19 09:57) [13]

ИМХО для решения уравнения методом деления отрезка пополам точность 1Е-30 не нужна, достаточно 1Е-6

← →
Dualsoft (2002-06-19 11:25) [14]

В данном случае желательно варьировать параметр
точности менять в зависимости от интевала поиска.
Отношение не интервала поиска к заданной точности
не должно превышать 1е19-1е20.

← →
Zavdim (2002-06-20 11:49) [15]

Вот я и говорю, что надо использовать относительную погрешность.

← →
Zavdim (2002-06-21 11:56) [16]

На самом деле надо бы выложить константы и полный текст проги.
Почему она все-таки циклилась?
Единственное существенное из вышесказанного это то, что нельзя дичиться с точностью, но в один момент abs(b-a) - должно было занулиться. Я поэсперементировал с нижеприведенным кодом, разница всегда уходила в ноль.

var a,b:extended;
procedure TForm1.FormCreate(Sender: TObject);
begin
a:=1e-30;
b:=1e+30;
l1.Caption:=FloatToStr(a);
l2.Caption:=FloatToStr(b);
l3.Caption:=FloatToStr(b-a)
end;

procedure TForm1.Button1Click(Sender: TObject);
var r:extended;
i:integer;
begin
r:=(a+b)/2;
i:=random(2);
if i=1 then a:=r else b:=r;
l1.Caption:=FloatToStr(a);
l2.Caption:=FloatToStr(b);
l3.Caption:=FloatToStr(b-a)
end;

Разобраться можно. Так что нечего дурить - давай свой код - посмотрим.

← →
Озадаченный (2002-06-21 12:45) [17]

Вот код программы:
procedure TForm1.Button1Click(Sender: TObject); var Nc,Nv,n,ni,Ec,Ev,Ed,Ea,EgT,Fa,Fx,x:Extended; r,a,b,k,Nd,dEd,Eg0,mdn,Na,dEa,alpha,mdp:Extended; T:integer; F:TextFile; begin k:=8.62e-5; Nd:=1e21; dEd:=0.044; Eg0:=1.21; mdn:=1.08; Na:=0; dEa:=0.044; alpha:=0.00028; mdp:=0.59; AssignFile(F,"otchet.txt"); Rewrite(F); for T:=1 to 1500 do begin EgT:=Eg0-alpha*T; Ec:=EgT/2; Ev:=-EgT/2; Ed:=Ec-dEd; Ea:=Ev+dEa; r:=T*T; // по другому почему-то нельзя возвести T в куб (T*T*T - не получается) r:=r*T; Nc:=4.831e21*sqrt(mdn*mdn*mdn)*sqrt(r); Nv:=4.831e21*sqrt(mdp*mdp*mdp)*sqrt(r); ni:=sqrt(Nc*Nv)*exp(-(EgT/(2*k*T))); a:=1e-30; //в принципе можно и от 1 использовать, вот только значения корней немного меняются почему-то b:=1e22; Fa:=a+Na/(2*exp((Ea-Ec-k*T*ln(a/Nc))/(k*T))+1)-ni*ni/a-Nd/(2*exp((Ec+k*T*ln(a/Nc)-Ed)/(k*T))+1); repeat x:=(a+b)/2; Fx:=x+Na/(2*exp((Ea-Ec-k*T*ln(x/Nc))/(k*T))+1)-ni*ni/x-Nd/(2*exp((Ec+k*T*ln(x/Nc)-Ed)/(k*T))+1); if Fa*Fx>0 then a:=x else b:=x; until (ABS(Fx)<1) or (ABS(b-a)<1e10); //ABS(b-a)<1 не получается посчитать n:=x; WriteLn(F,"n = "+FloatToStr(n)); end; CloseFile(F); end;

мне нужно найти n, но почему-то при больших T корень не ищется, хотя должен теоретически лежать в пределах от 1 до 1e21 максимум.

← →
Zavdim (2002-06-21 14:28) [18]

>r:=T*T; // по другому почему-то нельзя возвести T в куб (T*T*T - не получается)
>r:=r*T;

все просто, то что результат будет Extended он не знает заранее, а обрабатывает его как Integer, а теперь прикинь 1500^3 - это integer?

← →
Zavdim (2002-06-21 15:34) [19]

1. Ну тут штука такая, у меня в проге все к 0 сходилось и баста - хорошо, а здесь не сходилось - дествительно в последних разрядах прога висла и не меняла ни a ни b. Я ввел относительную погрешность - стала что-то делать. Но это ерунда.
2. А откуда ты уверен, что корни есть? У меня получилось что
на концах промежутка имеем одинаковые знаки. Ощущение такое, что
y=0 горизонтальная ассимптота. При b=1e22 получаются значения на пределе вычислительных возможностей Паскаля - это же чистый ноль.
3. Зачем константа Na=0, откуда Nd=1e21, при том, что к другим константам отношение более точное?
4. Возьми MathCad, Excel и построй графики этой функции.
5. Из цикла поиска мы выходим по двум причинам, в этой программе только по второй - область поиска съузилась, то что мы имеем не корни.
6. Я слегка модифицировал прогу, посмотри отчет и определись чего ты хочешь.

procedure TForm1.Button1Click(Sender: TObject);
var Nc,Nv,n,ni,Ec,Ev,Ed,Ea,EgT,Fa,Fb,Fx,x:Extended;
r,a,b,k,Nd,dEd,Eg0,mdn,Na,dEa,alpha,mdp:Extended;
T:integer;
TT:extended;
F:TextFile;
function ff(x:extended):extended;
begin
Result:=x+Na/(2*exp((Ea-Ec-k*TT*ln(x/Nc))/(k*TT))+1)-ni*ni/x-Nd/(2*exp((Ec+k*TT*ln(x/Nc)-Ed)/(k*TT))+1);
end;
begin
Button1.Enabled:=false;
k:=8.62e-5;
Nd:=1e21;
dEd:=0.044;
Eg0:=1.21;
mdn:=1.08;
Na:=0;
dEa:=0.044;
alpha:=0.00028;
mdp:=0.59;
AssignFile(F,"otchet.txt");
Rewrite(F);
for T:=1 to 1500 do begin
TT:=T;
EgT:=Eg0-alpha*TT;
Ec:=EgT/2;
Ev:=-EgT/2;
Ed:=Ec-dEd;
Ea:=Ev+dEa;
r:=exp(3*ln(TT));
Nc:=4.831e21*sqrt(mdn*mdn*mdn)*sqrt(r);
Nv:=4.831e21*sqrt(mdp*mdp*mdp)*sqrt(r);
ni:=sqrt(Nc*Nv)*exp(-(EgT/(2*k*TT)));
a:=1;
b:=1e+22;
Fa:=ff(a);
Fb:=ff(b);
repeat
x:=(a+b)/2;
Fx:=ff(x);
if Fa*Fx>0
then begin
a:=x;
Fa:=Fx;
end else
begin
b:=x;
end;
until (ABS(Fx)<1) or (ABS((b-a)/b)<1e-15);
n:=x;
WriteLn(F,"T="+IntToStr(T)+" <> n = "+FloatToStr(n)+" <> f(n)="+FloatToStr(n));
end;
CloseFile(F);
Button1.Enabled:=true;
end;

Желаю удачи. И вообще переходи на мыло, а то форум раздувается.

← →
Zavdim (2002-06-22 07:51) [20]

чтобы подвести некоторый итог:
1. Не надо пренебрегать точностью представления данных.
2. Прежде чем, что-то искать надо определиться с теорией, - ксожалению, по опыту знаю, что многие программисты плюют на математику и теорию, а это в итге порочно.
3. x^2-9=0 для приведенного приложения хороший пример, а вот, например, x^2-3x+2=0 - решений не даст, а даст огромные числа кои решениями не являются. Поэтому отделение корней не забава, а суровая необходимость.
4. Модификация должна состоять в том, что надо просто отделить корни (!!если они есть!!).

Трудно искать в темной комнате черную кошку, особенно...

Данные промахи полезны новичку, сам тонны соли на этом съел...
Но лучше самому разбираться в программе проходом по шагам, в данном случае ошибка выплывает при T=34.

Желаю удачи.

← →
Озадаченный (2002-06-24 11:59) [21]

Zavdim, Спасибо большое за помощь, попробую еще раз все перерасчитать и проверить точность

Вопрос Найти похожие ветки