Главная страница
    Top.Mail.Ru    Яндекс.Метрика
Текущий архив: 2004.10.24;
Скачать: CL | DM;

Вниз

Зацикливание программы   Найти похожие ветки 

 
saypn   (2004-10-08 22:30) [0]

Большая просьба, помогите разобраться с причиной зацикливания программы.

Есть задача: используя метод трапеции вычислить определенный интеграл от a до b для какой-либо функции.

Вывод формулы:
Надо найти площадь под графиком функции на отрезке [a,b], разбиваем его на n частей длины h=(b-a)/n, тогда площадь фигуры равна

S=(f(a)+f(a+h))/2*h+(f(a+h)+a(a+2h))/2*h+...+(f(a+(n-1)h)+f(a+nh))/2*h.

S=(f(a)+f(b))/2*h + сумма(i=1 до n-1) f(a+ih)*h

Вот таботающая программа:

program zad_pro_integral;

{$APPTYPE CONSOLE}

uses
 SysUtils;

type
 myfunc=function (x:double):double;

function integral(a,b:double;f:myfunc):double;
 var
   n,i:integer;
   eps,x,r0,r,h:double;
 begin
{n - количество трапеций, на которые разбивается фигура}
 n:=1000;
{r - вычесленная площадь фигуры}
 r:=0;
{eps - точность, с которой необходимо выполнить вычисление}
 eps:=1e-11;
 repeat
   r0:=r;
   n:=n*2;
   h:=(b-a)/n;
   r:=((f(b)+f(a)))/2*h;
   x:=a+h;
   for i:=1 to n-1 do begin
     r:=r+f(x)*h;
     x:=x+h;
   end;
 until abs(r-r0)<eps;
 result:=r;
 end;

function f1(x:double):double;
 begin
 result:=sin(x);
 end;

 var
s:integer;

begin
writeln(integral(0,pi/2,f1));
end.

Почему при значениях eps<1e-11 программа зацикливается. Чем это объясняется?

Огромная просьба откликнуться!


 
Algol   (2004-10-08 22:39) [1]


> Почему при значениях eps<1e-11 программа зацикливается.
> Чем это объясняется?

Это объясняется конечной точностью типа double, который ты используешь (15-16 значащих знаков).
Используй Extended, тогда точность повысится до 19-20 значащих знаков.


 
Юрий Зотов ©   (2004-10-08 22:55) [2]

Потому что такая точность оказывается недостижимой. А почему это происходит - здесь не расскажешь, это надо смотреть в книгах по вычислительной математике (точность метода, погрешность вычислений и ее накопление). Еще очень рекомендую прочитать вот эту статью:
http://www.delphikingdom.com/asp/viewitem.asp?catalogid=374

Простейшие варианты обхода - уменьшить шаг разбиения, либо использовать тип Extended (соответственно переключив FPU). Более сложные - либо шлифовать алгоритмы всех вычислений, чтобы повысить их точность (может и не помочь, все равно требуемая точность может быть и не достигнута), либо использовать методы более высокого порядка (что тоже не гарантирует такой высокой точности).

Лучший вариант - все вместе.


 
Algol   (2004-10-08 23:06) [3]


>  уменьшить шаг разбиения

Кажется вы невнимательно смотрели код. Шаг разбиения там и так уменьшается.


 
Юрий Зотов ©   (2004-10-08 23:29) [4]

> Algol   (08.10.04 23:06) [3]

Да, каюсь. Но на выводы это не влияет - из-за погрешностей вычислений алгоритм выходит на стационар и зацикливается.


 
Verg ©   (2004-10-08 23:36) [5]


> for i:=1 to n-1 do begin
>      r:=r+f(x)*h;
>      x:=x+h;
>    end;


Вот тут происходит накопление ошибки. Грубо говоря, суммарная ошибка прямо пропорциональна n, а оно (n), в свою очередь, растет в геометрической прогрессии (n := n*2) от итерации к итерации repeat until...


 
Юрий Зотов ©   (2004-10-08 23:49) [6]

> Verg ©   (08.10.04 23:36) [5]

Тут палка о двух концах. Как Вы и сказали, с увеличением n растет число проходов цикла интегрирования и накапливается ошибка, но с другой стороны уменьшается погрешность самого метода - а в итоге наступает такой момент, когда они компенсируют (или почти компенсируют) друг друга и алгоритм выходит на стационар (или квазистационар). Ситуация, для подобных программ типичная.


 
Verg ©   (2004-10-08 23:53) [7]


> [6] Юрий Зотов ©   (08.10.04 23:49)


Все верно. И чтобы я предложил в первую голову для "продления жизни до квазистационара", так это (кроме перехода на extended или смены метода интегрирования):

> for i:=1 to n-1 do begin
>      r:=r+f(a+h*i)*h;
>//////      x:=x+h;
>    end;
  x := x + h*n;


 
palva   (2004-10-09 00:13) [8]

Наверно, имелось ввиду
x := a+h*i
где i - номер шага.


 
Verg ©   (2004-10-09 00:34) [9]


> palva   (09.10.04 00:13)
> Наверно, имелось ввиду
> x := a+h*i


Да, хоть так.
Смысл один - избавится, где возможно, от накоплений в циклах. Погрешность вычисления f(x)*h накапливается в цикле + плюс сама погрешность вычисления х в этом же цикле вносит в интегральную сумму интегральную же погрешность. Так хоть от интегральной ошибки аргумента ф-ции избавиться - все вперед. :)



Страницы: 1 вся ветка

Текущий архив: 2004.10.24;
Скачать: CL | DM;

Наверх




Память: 0.47 MB
Время: 0.035 c
1-1097575940
Галинка
2004-10-12 14:12
2004.10.24
Как сделать не сортированный TStringList или TStrings


14-1096891606
Guest
2004-10-04 16:06
2004.10.24
Перестал пахать комп, пень первый 166 MHz


14-1097092784
saNat
2004-10-06 23:59
2004.10.24
Вопрос по C: как обратиться напрямую к биту


3-1095714619
Maxim______
2004-09-21 01:10
2004.10.24
тормоза BLOB в GDB


1-1097579697
DesWind
2004-10-12 15:14
2004.10.24
Аппроксимация





Afrikaans Albanian Arabic Armenian Azerbaijani Basque Belarusian Bulgarian Catalan Chinese (Simplified) Chinese (Traditional) Croatian Czech Danish Dutch English Estonian Filipino Finnish French
Galician Georgian German Greek Haitian Creole Hebrew Hindi Hungarian Icelandic Indonesian Irish Italian Japanese Korean Latvian Lithuanian Macedonian Malay Maltese Norwegian
Persian Polish Portuguese Romanian Russian Serbian Slovak Slovenian Spanish Swahili Swedish Thai Turkish Ukrainian Urdu Vietnamese Welsh Yiddish Bengali Bosnian
Cebuano Esperanto Gujarati Hausa Hmong Igbo Javanese Kannada Khmer Lao Latin Maori Marathi Mongolian Nepali Punjabi Somali Tamil Telugu Yoruba
Zulu
Английский Французский Немецкий Итальянский Португальский Русский Испанский