Ошибки в составном правиле Симпсона в C ++

Я пишу небольшую программу для аппроксимации интегралов на C ++ с использованием составного правила Симпсона. Кажется, это работает, но когда я вычисляю интеграл для разного количества интервалов (отсюда и цикл for), появляются некоторые ошибки, поскольку это кажется совершенно случайным. Вот мой код.

#include <iostream>
#include <cmath>

double f(double x)
{
    double y = x;
    return y;
}

int main()
{
    for(double n=1; n<20; n+=1)
    {
        double a = 1;
        double b = 4;
        double d = (b-a)/n;
        double i = a;
        double area = 0;
        while(i<b)
        {
            double j = i + d;
            area += (d/6) * ( f(i) + 4*f((i+j)/2) + f(j) );
            i = i + d;
        }
        std::cout << "Order " << n << " : " << area << std::endl;
    }
    return 0;
}

Результат выглядит следующим образом.

Order 1 : 7.5
Order 2 : 7.5
Order 3 : 7.5
Order 4 : 7.5
Order 5 : 7.5
Order 6 : 7.5
Order 7 : 9.30612
Order 8 : 7.5
Order 9 : 7.5
Order 10 : 8.745
Order 11 : 8.6281
Order 12 : 7.5
Order 13 : 7.5
Order 14 : 7.5
Order 15 : 7.5
Order 16 : 7.5
Order 17 : 8.22145
Order 18 : 8.18056
Order 19 : 7.5

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


c++ numerical-integration
person Naudin    schedule 24.01.2015    source источник
comment
вы не должны никогда использовать double в качестве индекса for цикла из-за ошибок округления. То же самое для while. В общем, вы должны проверять неравенство до некоторого конечного эпсилона, то есть std::abs(x - a) < eps. означает, что x достаточно близко к a, чтобы считаться равным.   -  person vsoftco    schedule 24.01.2015
comment
Я слишком ленив, чтобы самому подсчитывать, должен ли результат быть 7,5 в каждом случае?   -  person mmgross    schedule 24.01.2015

Ответы (1)


arrow_upward
1
arrow_downward

Вероятно, это вызвано сравнением чисел с плавающей запятой. 4.0<4.0 может быть правдой в этом случае для компьютера, из-за чего ваш цикл будет повторяться один раз слишком много.

Чтобы решить эту проблему, вы можете использовать epsilon для сравнения чисел:

#define EPS 1e-8
if(a<b-EPS){...} // this is analogous to < operator for ints
if(a<b+EPS){...} // this is analogous to <= operator for ints

Другой вариант - сгенерировать переменную i на лету, используя целые числа для управления циклом:

for(int k=0;k<n;k++){
    double i=a+(b-a)*k/n;
    // ...
}
person akrasuski1    schedule 24.01.2015