4

Пытаюсь вычислить частичную сумму знакопеременного ряда К/ (M*M), где К=1,-1,1,-1... Известно, что если М стремится к бесконечности, то сумма ряда равна (Пи*Пи)/12=0.8224670334241132... (точное значение невозможно, так как есть число Пи.) И даже получились очень симпатишные результаты!

Но почему-то точность падает на два порядка после 10 в -11 степени. По точности Пи и по разрядности переменных вроде есть запас.
В чем может быть моя ошибка - в типах переменных, в ключах вывода или в методе вычисления? Как это определить?

#include <windows.h>
#include <iostream>
#define _USE_MATH_DEFINES  // M_PI = 3.14159265358979323846
#include <math.h>
#include <iomanip>

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

3

3 ответа 3

3

В принципе, ничего неожиданного нет.

Вы вычисляете с типом 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, что, конечно, давало неправильный результат. Ликвидировав это проблему и подняв количество слагаемых, стало возможно добраться до такой точности, где и порядок слагаемых играет роль. Здесь вычисление от маленьких к большим наконец-то дало положительный эффект.

Ещё по теме: Алгоритм Кэхэна (здесь на русском).

10
  • Ох. Об этом я и не подумал. А правда, какой у вас максимальный размер intlong)? Если у вас int 32-разрядный, попробуйте вывести произведение 65000 * 65000, тогда станет понятна причина проблемы.
    – VladD
    Commented 26 окт. 2014 в 9:19
  • @qqqq1961: Я к тому, что у вас наверняка переполнение, и 65000 * 65000 «не вписывается» в int.
    – VladD
    Commented 26 окт. 2014 в 15:55
  • 1
    Уважаемые коллеги! Почему при суммировании в обратном порядке уже при ста лимонах членов получена ошибка НОЛЬ (гениальная догадка VladaD), а при прямом даже при 4-х миллиардах всего 10 в -14 и конца не видно?! Это мне напоминает, как меряли удава "от головы до хвоста..." И еще, ряд сходящийся из-за ошибки, что ли?! Или я нечаянно опроверг каку-нить аксиому параллельных, и они пересеклись? Куда бежать за нобелевкой???
    – qqqq1961
    Commented 26 окт. 2014 в 21:24
  • 1
    @qqqq1961: Сложение чисел типа double неассоциативно (хотя и коммутативно). Простой пример: если у нас действительные числа размером в 2 десятичных знака: (((1 + 0.07) + 0.06) + 0.05) + 0.04 = ((1 + 0.06) + 0.05) + 0.04 = (1 + 0.05) + 0.04 = 1 + 0.04 = 1 1 + (0.07 + (0.06 + (0.05 + 0.04))) = 1 + (0.07 + (0.06 + 0.09)) = 1 + (0.07 + 0.15) = 1 + 0.22 = 1.2 У вас происходит то же самое, только в бóльших масштабах. (У меня как раз описано в ответе.)
    – VladD
    Commented 26 окт. 2014 в 21:38
  • 1
    @qqqq1961: При прямом суммировании, думаю, начиная с некоторого слагаемого при прибавлении сумма вообще не меняется.
    – VladD
    Commented 26 окт. 2014 в 21:46
2

Суммировать следует с наименьших чисел, чтобы свести к минимуму накапливаемую погрешность. Посмотри мою небольшую лекцию.

Следующая задачка: подсчитать сумму 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) это, наверное, неактуально.

2
  • Вот теперь всё понятно!!! Спасибо, Валентин, за потраченное на меня время, наверняка и другим будет полезно такое разъяснение. Осталось только понять главное. Вот Вы говорите: "...так будет аккуратнее. Но для 1/(m*m) это, наверное, неактуально." Как же узнать - когда и как аккуратнее и актуальнее?! Как сравнивать разные алгоритмы и с чем. Неужели IEEE-754 даёт гарантии для бесконечных рядов и иррациональных чисел?
    – qqqq1961
    Commented 28 окт. 2014 в 18:25
  • Суммировать разности стоит для очень близко стоящих друг к другу соседних элементов. Просто (эмпирически так предполагаю) не стоит тратить дополнительные усилия на усложнение кода для данной задачи. IEEE-754 не дает гарантий абсолютной точности. Можно лишь постараться свести к минимуму издержки дискретной арифметики
    – Valentin
    Commented 29 окт. 2014 в 1:16
1

@qqqq1961, боюсь тут все сложнее, чем Вам кажется.

Внимательно посмотрите на простую программу, ее результаты и сравненние их с подсчетом в bc с точностью 40 знаков.

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h>

#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 на самом-то деле ограничена точностью аппаратуры.

Отсутсвие ПРАВИЛЬНОГО ОКРУГЛЕНИЯ и приводит к странным результатам (особенно при подсчете суммы, начиная с больших членов).

8
  • Да, заметил. >Заметили, что доверять можно только 16 десятичным знакам Думаю, что столько знаков после умножения не нужно, если исходные 16 знаков и сумма ряда, видимо, тоже не должна превышать 16. Даже если сопроцессор, как я читал, использует 80. >3.1415926535897931 * 3.1415926535897931 = 9.86960440108935774884804450080761
    – qqqq1961
    Commented 27 окт. 2014 в 12:53
  • @qqqq1961, да это все к вопросу -- насколько можно доверять* Limit = (M_PI*M_PI)/12; ... Summa = ... ... Error = Limit - Summa; вот таким вычислениям? Понятно же, что младшие разряды как у Limit, так и у Summa недостоверны.
    – avp
    Commented 27 окт. 2014 в 13:09
  • Понятно же, что младшие разряды как у Limit, так и у Summa недостоверны. Эй, эй! Что значит недостоверны?! После какой степени ошибки данные вычисления врут? 10 в-4? 10 в-14? Вы говорите вообще про точность вычисления суммы, даже при, скажем, 1000 членов или только про близкие к концу double значения? И что конкретно плохого в: Limit = (M_PI*M_PI)/12; Summa = ... Error = Limit - Summa;? Как еще можно считать если не нужна предельная дочность?
    – qqqq1961
    Commented 27 окт. 2014 в 18:07
  • @qqqq1961, к сожалению я не могу конкретно ответить на Ваши вопросы. -- Кое-что можно прочесть здесь (или нагуглить другие материалы...). Также можно надеяться, что кто-нибудь на ХК достаточно компетентен и не ленив... В любом случае, позвольте привести еще один "отрицательный" пример: double lim1 = (M_PI * M_PI) / 12, lim2 = (M_PI / 12) * M_PI; printf ("lim1 = %.16f lim2 = %.16f\n", lim1, lim2); lim1 = 0.8224670334241132 lim2 = 0.8224670334241131 Ну и с каким из lim сравнивать правильно? (подозреваю, что со вторым).
    – avp
    Commented 27 окт. 2014 в 21:48
  • @qqqq1961, @avp: Угу, и вот ещё туда же классическая статья (на неё есть ссылка в статье на Хабре).
    – VladD
    Commented 27 окт. 2014 в 23:04

Ваш ответ

Нажимая «Отправить ответ», вы соглашаетесь с условиями пользования и подтверждаете, что прочитали политику конфиденциальности.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.