复合辛普森规则代码不断得到错误的输出？

Question

算法：enter image description here

我正在尝试执行此复合辛普森规则，该规则将计算以下积分：1/sqrt(x)，其结果应为：2

但是，我一直收到错误的输出，例如 1.61663

#include <cmath>
#include <iostream>
using namespace std;

double f(double n)
{

    return 1/sqrt(n);
}

double simpson(double a, double b, double n)
{

double x0=f(a)+f(b);
double h=(b-a)/(n);
double x1=0,x2=0;
double x=0;
for(int i = 1 ; i <n;i++){
        x=a+(i*h);
        if(i%2==0)
        {
            x2=x2+f(x);
        }
        else
        {
            x1=x1+f(x);
        }
    }
        x1=(h*(x0+2*x2+4*x1))/3;
        return x1;
}

 int main(){

            cout<<"Integral is: "<<" "<<simpson(0.0004,1,20)<<" "<<endl;
    }

Answer 1

问题不在您的代码中，而在您尝试集成的函数中。随着 x 趋于零，此函数发散到无穷大。衍生品也有分歧。

对于 a 较小的区间 [a, 1]，误差项以 O[1/(N^4 a^4.5)] 为界。这就是为什么要计算这个区间内的积分，网格应该非常密集以获得合理的误差范围。

simpson(0.0004, 1, N) 生成以下值：

N      Result         Error
--------------------------------
20     2.549041009    0.5890
200    1.986462457    0.0265
2000   1.960181808    1.8181e-4
20000  1.960000049    4.9374e-8 
200000 1.960000000    5.0810e-12

确实，对于大 N，我们越来越接近精确值 1.96，但存在错误 O(1/N^4)。

复合辛普森规则代码不断得到错误的输出？

Composite simpson rule code keep getting wrong output?

c++

algorithm

math

calculus