#visual_cpp
Пытаюсь вычислить частичную сумму знакопеременного ряда К/ (M*M), где К=1,-1,1,-1... Известно, что если М стремится к бесконечности, то сумма ряда равна (Пи*Пи)/12=0.8224670334241132... (точное значение невозможно, так как есть число Пи.) И даже получились очень симпатишные результаты! Но почему-то точность падает на два порядка после 10 в -11 степени. По точности Пи и по разрядности переменных вроде есть запас. В чем может быть моя ошибка - в типах переменных, в ключах вывода или в методе вычисления? Как это определить? #include#include #define _USE_MATH_DEFINES // M_PI = 3.14159265358979323846 #include #include int main() { int N=1, M; double Summa, Error, Limit, One; cout << "Вычисление частичной суммы ряда и сравнение её с пределом.\n\n" << endl << endl; Limit = (M_PI*M_PI)/12; while (N < 65000) { Summa=0; Error=0; One=1; for ( M = 1; M < N; M++) { Summa += ( One / (M*M) ); // здесь вычисляем сумму первых членов One=-One; // ряд знакопеременный } Error = Limit - Summa; cout << "Cумма " << N << " членов ряда = " << setprecision(10) << Summa << "\t Ошибка = " << setprecision(3) << abs(Error) << endl; N=N+1000; // шаг изменения кол-ва членов частичной суммы } cout << " Предел суммы ряда: " << setprecision(16) << Limit << endl; } Cумма 1001 членов ряда = 0.8224665339 Ошибка = 5e-007 - больше членов - меньше ошибка Cумма 2001 членов ряда = 0.8224669085 Ошибка = 1.25e-007 Cумма 3001 членов ряда = 0.8224669779 Ошибка = 5.55e-008 Cумма 4001 членов ряда = 0.8224670022 Ошибка = 3.12e-008 Cумма 5001 членов ряда = 0.8224670134 Ошибка = 2e-008 Cумма 6001 членов ряда = 0.8224670195 Ошибка = 1.39e-008 Cумма 7001 членов ряда = 0.8224670232 Ошибка = 1.02e-008 Cумма 8001 членов ряда = 0.8224670256 Ошибка = 7.81e-009 Cумма 9001 членов ряда = 0.8224670273 Ошибка = 6.17e-009 Cумма 10001 членов ряда = 0.8224670284 Ошибка = 5e-009 ... Cумма 53001 членов ряда = 0.8224670333 Ошибка = 1.29e-010 Cумма 54001 членов ряда = 0.8224670333 Ошибка = 1.03e-010 Cумма 55001 членов ряда = 0.8224670334 Ошибка = 7.2e-011 Cумма 56001 членов ряда = 0.8224670334 Ошибка = 3.42e-011 Cумма 57001 членов ряда = 0.8224670334 Ошибка = 1.24e-011 Cумма 58001 членов ряда = 0.8224670335 Ошибка = 7.14e-011 Cумма 59001 членов ряда = 0.8224670336 Ошибка = 1.49e-010 <- и вдруг... Cумма 60001 членов ряда = 0.8224670337 Ошибка = 2.54e-010 Cумма 61001 членов ряда = 0.8224670338 Ошибка = 4.06e-010 Cумма 62001 членов ряда = 0.8224670341 Ошибка = 6.43e-010 Cумма 63001 членов ряда = 0.8224670345 Ошибка = 1.07e-009 <- падает точность... Cумма 64001 членов ряда = 0.8224670355 Ошибка = 2.05e-009 Предел суммы ряда: 0.8224670334241132
Ответы
Ответ 1
В принципе, ничего неожиданного нет. Вы вычисляете с типом double, который имеет около 16 значащих цифр. Теперь, если само слагаемое имеет порядок 1e-5 * 1e-5 = 1e-10, а вы складываете его с частичной суммой, которая имеет величину порядка единицы, дальнейшие разряды теряются, и реально прибавляются лишь последние 6 десятичных цифр: част. сумма: x.xxxxxxxxxxxxxxxx слагаемое: 0.0000000000xxxxxxxxxxxxxxxx ------------------------------------------ результат: x.xxxxxxxxxxxxxxxx То есть, точность слагаемого при таком вычислении фактически ограничена 1e-16. Когда у вас количество слагаемых порядка 1e+5, суммарная точность получается порядка 1e-11. Попробуйте вычислять с конца: сначала наименьшие слагаемые. Тогда частичная сумма будет, возможно, не так сильно понижать точность. (Не уверен, что это поможет.) Обновление: на самом деле тут были две проблемы. Первая — переполнение int: произведение M * M не помещалось в 32-битный int, поэтому вычислялось по модулю 2^32, что, конечно, давало неправильный результат. Ликвидировав это проблему и подняв количество слагаемых, стало возможно добраться до такой точности, где и порядок слагаемых играет роль. Здесь вычисление от маленьких к большим наконец-то дало положительный эффект. Ещё по теме: Алгоритм Кэхэна (здесь на русском).Ответ 2
@qqqq1961, боюсь тут все сложнее, чем Вам кажется. Внимательно посмотрите на простую программу, ее результаты и сравненние их с подсчетом в bc с точностью 40 знаков. #include#include #include #define __S_PI(x) # x #define _S_PI(x) __S_PI(x) #define S_PI _S_PI(M_PI) int main(int ac, char *av[]) { printf ("M_PI = %.16f M_PI ** 2 = %.16f (math.h: %s)\n", M_PI, M_PI * M_PI, S_PI); return puts("End") == EOF; } avp@avp-xub11:hashcode$ gcc c.c avp@avp-xub11:hashcode$ ./a.out M_PI = 3.1415926535897931 M_PI ** 2 = 9.8696044010893580 (math.h: 3.14159265358979323846) End avp@avp-xub11:hashcode$ bc bc 1.06.95 Copyright 1991-1994, 1997, 1998, 2000, 2004, 2006 Free Software Foundation, Inc. This is free software with ABSOLUTELY NO WARRANTY. For details type `warranty'. scale=40 3.1415926535897931 * 3.1415926535897931 9.86960440108935774884804450080761 3.14159265358979323846 * 3.14159265358979323846 9.8696044010893586188178821328931344231716 avp@avp-xub11:hashcode$ Заметили, что доверять можно только 16 десятичным знакам (15 после запятой)? Особо посмотрите, как падает точность после умножения (M_PI * M_PI). Другими словами, точность вычислений с double на самом-то деле ограничена точностью аппаратуры. Отсутсвие ПРАВИЛЬНОГО ОКРУГЛЕНИЯ и приводит к странным результатам (особенно при подсчете суммы, начиная с больших членов). Ответ 3
Суммировать следует с наименьших чисел, чтобы свести к минимуму накапливаемую погрешность. Посмотри мою небольшую лекцию. Следующая задачка: подсчитать сумму 10 000 000 элементов ряда: S = 1 + 1/2 + 1/3 + 1/4 + 1/5 + 1/6 + . . . + 1/10000000. По идее от перемены слагаемых сумма не меняется, т.е. если буду складывать в обратном порядке: S = 1/10000000 + 1/9999999 + 1/9999998 + . . . + 1/3 + 1/2 + 1 то сумма должна быть такой же. Это верно для классической арифметики. Но мы имеем дело с компами, в которых заложена дискретная арифметика. Напишем нижеследующий код и скомпилируем его: program Summa; var s:single; i:longint; begin s:=0.0; for i:=1 to 10000000 do s:=s+1.0/i; //подсчет суммы S = 1 + 1/2 +...+ 1/10000000 Writeln(s:0:2); //вывод результата на экран s:=0.0; for i:=10000000 downto 1 do s:=s+1.0/i; //подсчет суммы S = 1/10000000 +...+ 1/2 + 1 Writeln(s:0:2) //вывод результата на экран end. На экране будут такие результаты: 15.40 16.69 Результаты разные. В чем дело? Какая сумма из них более достоверна? И какая отсюда мораль? Проблема заключается в конечности представления чисел в компе. Для данных типа single отводится 4 байта, в которых 1 бит отведен под знак, 8 бит – под двоичную степень, а оставшиеся 23 бита – под мантиссу ==> точность представления чисел ограничена 23 битами или 7-8 значащими десятичными цифрами. Т.е. в комп можно еще записать число 7.7777777 или 7.7777778, но никак не 7.7777777000000000007777777 – такое число просто округляется/обрезается до 23 бит для мантиссы. Для того, чтобы проще объяснить работу нашего вышеприведенного примера, допустим, что точность чисел ограничена 2 значащими числами. Тогда в цикле подсчета суммы в прямом направлении имеем: Проход 1: s := 0.0 + 1.0/1.0 = 1.0 Проход 2: s := 1.0 + 1.0/2.0 = 1.5 Проход 3: s := 1.5 + 1.0/3.0 = 1.5 + 0.33333333… --> 1.8 (округление до 2 цифр) . . . Проход 100: s := s + 1.0/100.0 = s + 0.01… --> s (округление и потеря младших цифр) Проход 101: s := s + 1.0/101.0 = s + 0.0099… --> s (округление и потеря младших цифр) . . . При такой точности после 100-го прохода (фактически еще раньше) сумма в цикле будет адекватной s := s + 0.0 из-за выбрасывания младших цифр в накапливаемой сумме. Если же проделывать в обратном направлении: Проход 1: s := 0.0 + 1.0/10000000.0 = 0.0000001 (результат сохраняется до 2 значащих цифр) Проход 2: s := 0.0000001 + 1.0/9999999.0 = 0.0000001+ 0.00000010000001…--> 0.0000002 Проход 3: s := 0.0000002 + 1.0/9999998.0 = 0.0000001+ 0.00000010000002…--> 0.0000003 . . . то мы получим более достоверную сумму. Кстати, если декларировать для s тип уже не single (одинарной точности), а double (двойной точности – 8 байт), то на экран будет выведено следующее: 16.70 16.70. Но если же суммировать для бOльшего числа элементов (к примеру, до 100000000000000 элементов), то опять могут возникнуть искажения. Мораль: для суммирования данных с большим разбросом значений подсчет суммы следует начинать с наименьших элементов, для того чтобы свести к минимуму потерю информации из-за дискретной природы данных, представляемых в компе. И еще что бы посоветовал - раз ряд знакопеременный, то обычно стоило бы суммировать разности соседних элементов, так будет аккуратнее. Но для 1/(m*m) это, наверное, неактуально.
Комментариев нет:
Отправить комментарий