Ошибка подсчета интеграла при ОЧЕНЬ маленьком шаге

Ниже приведен алгоритм вычисления интеграла методом прямоугольника:



#include <iostream.h>

#include <iomanip.h>



double Parabola(double x) { return x*x; }



int main(void)

{

double a = 0.0, b = ???, dx = ???, square = 0.0;



for(double x = a; x <= b; x += dx)

square += dx*Parabola(x);

cout<<setprecision(20)<<square<<endl;



return 0;

}



Имеются следующие результаты:



a = 0.0; b = 1000.0 dx = 0.001;

square = 333333833.3244913816



a = 0.0; b = 1000.0 dx = 0.0001;

square = 333333383.2359925508



a = 0.0; b = 1000.0 dx = 0.00001;

square = 333333337.3506906628



a = 0.0; b = 1000.0 dx = 0.000001;

square = 333333334.242036581



a = 0.0; b = 1000.0 dx = 0.0000001;

square = 333333423.9663200378



a = 0.0; b = 1000.0 dx = 0.00000005;

square = 3333331??.?????????? // точно не помню







a = 0.0; b = 10.0 dx = 0.001;

square = 333.3833350000317637



a = 0.0; b = 10.0 dx = 0.0001;

square = 333.3383333496889804



a = 0.0; b = 10.0 dx = 0.00001;

square = 333.3338333242379008



a = 0.0; b = 10.0 dx = 0.000001;

square = 333.3333833585027719



a = 0.0; b = 10.0 dx = 0.0000001;

square = 333.3333385448191848



a = 0.0; b = 10.0 dx = 0.00000001;

square = 333.3333209307022571



a = 0.0; b = 10.0 dx = 0.000000001;

square = 333.3333057475213082



А вот вопрос. Почему при уменьшении шага dx точность начинает ухудшаться?

 

Ale

Это обычное явление при применении численных методов. У нас есть два источника неточностей: 1) ошибки самого алгоритма и 2) ошибки в данных и их последствия(округление и т.д.). С уменьшением шага 1) уменьшается а 2) возрастает и начиная с некоторого момента перевешивает. Оптимальный шаг получается как правило в районе О(sqrt(е)), где е - машинная точность, обычно 10^-16. Что твоими результатами лишний раз и подтверждается.



Таков общий принцип, при желании можно наверное и точное выражение для ошибки вывести. Для численного дифференцирования мы, помнится, на какой-то контрольной выводили, ну а интегрирование вроде не сильно отличается..

 

Я не помню метода прямоугольника для численого интегрирования, может это метод трапеций? А вообще лучше оперировать алгоритмами в виде блок-схем, назвать название метода.

 

То, что точность ограничена это понятно. Поэтому ес-но нужно работать с 80-битными double extended и не делать умножения на константу dx в цикле - это лишь ухудшает точность. Но по всей видимости тут замешиваются и побочные эффекты. Нельзя однозначно утверждать, что в приведенном диапазоне значений точность ухудшается, т.к. метод то сам по себе неточный: что значит x <= b ? Это значит что мы всегда завышаем результат и при определенных dx это завышение компенсирует отрицательную неточность вычисления. Если к примеру взять условие x < b, то при 80-битных вычислениях точность в указанном диапазоне dx только увеличивается. Особенно наглядно случайная компенсация погрешности проявляется в методе трапеций (оличается от приведенного варианта только тем, что при x=b берется x*x/2). В этом случае при dx = 0.000001 наблюдается якобы "резкое увеличение точности" на 5-6 порядков, а затем все возвращается назад и точность монотонно растет с уменьшением dx. Другой эффект наблюдается при аппроксимации прямоугольниками по средней точке интервала. Тут точность сначала растет, затем падает. Но это лишь следствие того, что при определенных dx в среднем лучше компенсируются ошибки, обусловленные +\- отклонениями параболы от апроксимирующего прямоугольника.

Код (Text):

1. Прямоугольник (for x=a,x<b)
2. 333.33[b]8[/b]333349999788, dx = 0.0001
3. 333.333[b]8[/b]33333499990, dx = 0.00001
4. 333.3332[b]8[/b]3333361439, dx = 0.000001
5. 333.33333[b]8[/b]333088254, dx = 0.0000001
6. 333.333333[b]8[/b]30981254, dx = 0.00000001
7. 333.3333332[b]5[/b]8103427, dx = 0.000000001
8.  
9. Трапеция (for x=0,x<b плюс b*b/2)
10. 333.3[b]4[/b]3333450000288, dx = 0.0001
11. 333.33[b]4[/b]333334499991, dx = 0.00001
12. 333.3333333333[b]6[/b]1439, dx = 0.000001  <- случайное совпадение
13. 333.3333[b]4[/b]3333088354, dx = 0.0000001
14. 333.33333[b]4[/b]330981255, dx = 0.00000001
15. 333.3333333[b]0[/b]8103427, dx = 0.000000001
16.  
17. Прямоугольник по средней точке (for x=a+dx/2,x<b)
18. 333.3333333[b]2[/b]4999789, dx = 0.0001
19. 333.333333333[b]2[/b]49992, dx = 0.00001
20. 333.3333333333[b]5[/b]8937, dx = 0.000001
21. 333.333333333[b]0[/b]88250, dx = 0.0000001
22. 333.33333333[b]0[/b]981342, dx = 0.00000001
23. 333.3333333[b]0[/b]8104280, dx = 0.000000001

 

По большому счету dx должно быть разным для каждого узла и зависить от величины функции в этом узле.

 

leo

Никак не могу понять, по какому принципу ты цифры жирным выделял, почему например в

выделена 8 а не 2, в

5 а не 2?

Да, таким образом мы немного передвинем проблему и критичной станет область О(10^-10). А с прямоугольниками по средней точке интересно получается, ты не смотрел что будет дальше, для еще меньших dx?

 

Для повышения точности можно заменить суммирование с dx на умножение:

Код (Text):

1.  
2. double Sum=0;
3. __int64 n=(b-a)/dx,i;
4. for (i=0; i<=n; i++)
5.  Sum += y(a+i*dx);
6.  

Тогда получаются такие результаты (a=0, b=10, второй столбик - ошибка, третий - dx):



333.32833335000032090000 4.999983e-03 1.000000e-04

333.33283333349100980000 4.999998e-04 1.000000e-05

333.33338333334853590000 5.000002e-05 1.000000e-06

333.33333833338593880000 5.000053e-06 1.000000e-07

333.33333283319836940000 5.001350e-07 1.000000e-08

333.33333327232583090000 6.100749e-08 1.000000e-09



Не знаю, что там дальше будет, но пока ошибка только падает.



_DEN_

Это как?

 

Stiver



> "Никак не могу понять, по какому принципу ты цифры жирным выделял"



Вообще-то на шару, но можно сказать - по порядку ошибки

При значении ..332833.. ошибка будет ..000500..



> "Да, таким образом мы немного передвинем проблему и критичной станет область О(10-10)"



Ес-но. Причину ошибки уже Artem назвал - это накопление ошибки при x=x+dx. До меня это сразу дошло, как только проветрил мозги никотином, да вот связь со Златоглавой оборвалась по известным причинам

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

А если делать умножение или для ускорения использовать комбинированный метод (например 1000..10000 складываем, затем берем уточненное значение умножением и т.д.) то ес-но получается точнее. Вот пример с 80-битными числами, метод трапеций с суммированием x по 1000 dx, затем уточнение x путем умножения:

Код (Text):

1. 333.3333333[b]5[/b]0000010, dx = 0.000100000
2. 333.333333333[b]4[/b]99936, dx = 0.000010000
3. 333.33333333333[b]5[/b]049, dx = 0.000001000
4. 333.3333333333333[b]5[/b]6, dx = 0.000000100
5. 333.333333333333[b]2[/b]15, dx = 0.000000010
6. 333.33333333333[b]4[/b]080, dx = 0.000000001

Вот только не проверил - можно ли уменьшить погрешность, если взять последний dx'=|x-b|. Еще интересно, не лучше ли будет "уточнить" значение dx, как dx"=(b-a)/n где n =round((b-a)/dx) или это вещь в себе ?

 

Резюме:

Если аргумент функции y(x) считается методом приращения x=x+dx, то при малых dx и соответственно большом числе слагаемых N=(b-a)/dx из-за ограниченной точности двоичного представления dx происходит нарастание погрешности расчета x. Соответственно, по условию (x < b) или (x <= b) расчет останавливается на значении x не равном b-dx или b и мы теряем последний прямоугольник (трапецию) шириной b-x.

Соответственно даже при достаточно больших dx максимальная ошибка может достигать y(b)*dx если конечное значение x чуть-чуть перескакивает за b. И это не погрешность метода прямоугольников или трапеций, а влияние ошибок округленя dx при суммировании x. Но могут быть и "удачные" значения dx при которых конечное значение x < b и b-x << dx.

Ну а при очень малых dx и соответсвенно больших N проявляются ошибки округления при суммировании s. Величина ошибки получается порядка y_ср*N*e/2 = s*e/(2dx), где e=1/(2^m), m - число бит мантиссы (54 для double и 64 для extended). Для dx=1E(-9) и s=333 получаем ошибку порядка 1E(-5) для double и 1E(-8) для extended, что примерно и наблюдаем.



Ес-но напрашивается вариант уточнения результата путем добавления недостающего интервала. Т.е. производим обычное суммирование с интервалом dx до x <= b-dx, а затем учитываем последний интервал dx'=b-x.

Код (Text):

1. Обычные трапеции, for x=0,x<b,x+=dx {s=s+x*x}; s=(s+b*b/2)*dx;
2. 333.3[b]4[/b]3333349999788, dx = 0.0001  <- завышение т.к. последний интервал dx' < dx
3. 333.33[b]4[/b]333333499990, dx = 0.00001
4. 333.3333333333[b]6[/b]1439, dx = 0.000001 <- удачное совпадение
5. 333.3333[b]4[/b]3333088254, dx = 0.0000001
6. 333.33333[b]4[/b]330981254, dx = 0.00000001
7. 333.3333333[b]0[/b]8103427, dx = 0.000000001 <- ошибки накопления суммы
8. Трапеции с уточнением последнего интервала
9. for x=0,x<=b-dx,x+=dx {s=s+x*x}; s=(s+x*x/2)*dx+(x*x+b*b)*(b-x)/2;
10. 333.3333333[b]4[/b]9999819, dx = 0.0001  <- ошибка метода уменьшается с уменьшением dx
11. 333.333333333[b]5[/b]06176, dx = 0.00001
12. 333.3333333332[b]9[/b]5099, dx = 0.000001
13. 333.333333333[b]1[/b]18328, dx = 0.0000001  <- ошибка округления суммы растет с уменьшением dx
14. 333.33333333[b]5[/b]285813, dx = 0.00000001
15. 333.3333332[b]9[/b]3391519, dx = 0.000000001

Результаты, как говорится, на лицЕ

 

leo

Хм, не совсем понимаю, что ты имеешь в виду. На всякий случай подведу итоги темы, а заодно и некую теоретическую базу под наши рассуждения



Опытным путем для всех трех алгоритмов было установлено, что точность результата сначала монотонно растет, а начиная с некоторого момента так же монотонно падает. Максимальная точность достигается в случае прямоугольника - при h=10^-8-10^-9, в случае трапеции и прямоугольника по средней точке - при h=10^-5 (все числа приблизительные, h:=dx).



Вопрос: почему такое происходит и как если дан алгоритм, найти оптимальное h ?



Как уже указывалось, нас интересуют два вида ошибок. Начнем с алгоритмических("ошибка дискретизации") А:

алгоритм прямоугольника - А=О(h)

алгоритмы трапеции и прямоугольника п.с.т. - А=О(h²) (отличаются только константой)

Доказательства приводить не буду т.к. тривиальны, да и в сети найти можно.



Перейдем теперь к ошибкам Е, связанным с конечной точностью компьютера. Все три алгоритма сводятся к вычислению некоторой суммы. Каждое слагаемое известно нам с точностью до е ("машинная точность") - O(e). Будем считать, что мы работаем с double, а значит е=10^-16. Количество слагаемых равно О(1/h), отсюда Е=О(е/h).



Суммарная ошибка F равна естественно Е + А. Все факторы, константы и прочее выкидываем для наглядности.

1) алгоритм прямоугольника:

F = O(h) + O(e/h) = h + e/h

Это выражение становится минимальным при h=sqrt(e)=10^-8, ошибка тогда равна тоже sqrt(e)=10^-8.



2) алгоритмы трапеции и прямоугольника п.с.т.:

F = O(h²) + O(e/h) = h² + e/h

Это выражение становится минимальным при h=e^(1/3)=10^-5, ошибка тогда равна e^(2/3)=10^-10.



Как видим, теория полностью согласуется с нашими наблюдениями и наоборот. Что интересно, попутно мы вычислили и максимально достижимую точность для каждого алгоритма(при использовании double конечно). Т.е. например с помощью прямоугольника больше чем 7-8 значащих цифр получить не удасться, что бы мы не делали.



Вся эта теория относится естественно только к самим алгоритмам в чистом виде. Понятно что кривая реализация может добавить сколько угодно дополнительных ошибок, пример: условие x < b.



В дополнение и мои результаты(a=0, b=1):.



Прямоугольник:

0.33283350000000017000 : 0.001

0.33328333500003693000 : 0.0001

0.33332833335128759000 : 0.00001

0.33333283332608060000 : 0.000001

0.33333328348568303000 : 0.0000001

0.33333332688968065000 : 0.00000001

0.33333334009230658000 : 0.000000001

0.33333330578792780000 : 0.0000000001



Трапеция:

0.33333349999999956000 : 0.001

0.33333333500003698000 : 0.0001

0.33333333335128712000 : 0.00001

0.33333333332608051000 : 0.000001

0.33333333348565591000 : 0.0000001

0.33333333188961028000 : 0.00000001

0.33333334059217068000 : 0.000000001

0.33333330583793208000 : 0.0000000001



Прямоугольник п.с.т.:

0.33333324999999953000 : 0.001

0.33333333250003599000 : 0.0001

0.33333333332628828000 : 0.00001

0.33333333332583109000 : 0.000001

0.33333333348565325000 : 0.0000001

0.33333333188961017000 : 0.00000001

0.33333334059217068000 : 0.000000001

0.33333330583793208000 : 0.0000000001





Artem

От абсолютной ошибки толку как правило мало, в таких случаях надо давать или относительную или количество значащих цифр(что в принципе одно и то же)



Добавил: писал, когда последнего сообщения leo еще не было. При суммировании учитывал последний интервал h'=b-x.

 

Artem

А вот так. Операции + и - чувствительны к разности порядков операндов. Поэтому точность и зависит от значений и поведения функции.

 

>Stiver

Помнишь о Обобщенной Функции Миттаг-Лефлера нечто подобное происходило и там

 

Попробуйте прооптимизировать этим.

Ака найти золотую середину.



_62600682__OCHKP.7z



Там на листе Y в Xmin Xmax пишешь задаваемые параметры (4 штуки)

А в таблице X1 X2 X3 X4 Y в Y подставляешь высчитанное по показываемым X1 X2 X3 X4 значение ошибки твоей функции.

 

murtix

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

 

Stiver

Вот если бы с возрастанием аргумента количество слагаемых тоже увеличивалось(экспоненциально),



А ведь увеличивается, чем больше аргумент, тем больше слагаемых необходимо взять



И по ходу вопрос:

Имеется несобственный интеграл от 0 до +бескон. какое должно быть условие окончания вычислений.

 

murtix

Хм.. насколько мне известно, никто несобственные интегралы численно "в лоб" не вычисляет. Так как понятно, что такого условия в общем виде нет и быть не может. Обычно переменную трансформируют сначала таким образом, что интервал интеграции становится конечным и потом работают уже с обычными методами. Системы типа Maple пытаются правда, пользуясь возможностью производить символьные вычисления, этот процесс автоматизировать, но пока с очень переменным успехом

 

Stiver

Но ведь бескон. (сходящиеся, конечно) ряды вычисляют, с точки зрения вычисления разница по моему не очень большая, я делал так: брал отрезки (0,1)+(1,2)+(2,3)+....

и останавливал процесс когда очередной член становился достаточно малым по отношению к сумме. Вроде все было достаточно правдоподобно. Разве это неправильно?

 

murtix

А как определяют, что ряд сходящийся? Скажем на примере ряда 1/n ?

 

Ну, это уже математика, есть же условия схождения. (Коши-Доламбера и т.п.)

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

и насколько я помню ряд 1/(n^p) сходится в случае р>1 и не сходится при р<=1.

 

murtix

Погоди, не вали все в одну кучу Если я правильно понял, речь шла о _численном_ интегрировании. Т.е. задача в обшем виде выглядит так: дана функция f (причем дана в виде черного ящика!) и границы a и b, между которыми ее надо интегрировать. При такой постановке вопроса определить, сходяшийся интеграл или нет, невозможно. Соответственно нельзя определить и когда заканчивать вычисления.



Другое дело когда функцию можно изучить с помощью "математики". Если интеграл сходящийся, тогда действительно существует такое w, что для x>w можно применить твой метод

Но проблема состоит в том, что w опять же невозможно вычислить численно, только через теорию(зависит помимо функции и от требуемой точности результата). Если брать просто все время w=0 как у тебя, то ничего хорошего не выйдет(если конечно случайно не попадешь на функцию у которой действительно w=0 )



Все это относится естественно и к бесконечным рядам.