Зацикливание программы

← →
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 растет число проходов цикла интегрирования и накапливается ошибка, но с другой стороны уменьшается погрешность самого метода - а в итоге наступает такой момент, когда они компенсируют (или почти компенсируют) друг друга и алгоритм выходит на стационар (или квазистационар). Ситуация, для подобных программ типичная.

> [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 - номер шага.

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

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

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