численное интегрирование — ошибки в правиле составного Симпсона в переполнении стека
Я пишу небольшую программу для аппроксимации интегралов в 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
Я думаю, что это может быть связано с типами переменных или сохранением проблем, хотя я не могу понять, в чем проблема. Любая помощь приветствуется.
Решение
Это, вероятно, вызвано сравнением чисел с плавающей запятой. 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;
// ...
}
Другие решения