复合辛普森规则无穷大输出

Posted 2023-02-22

【中文标题】复合辛普森规则无穷大输出

【英文标题】：Composite simpson rule infinity output

【发布时间】：2019-12-08 16:19:39

【问题描述】：

算法：

我有以下代码使用复合辛普森规则计算不正确的积分，我正在尝试评估积分exp(-x)/sqrt(1-x)，其中a= 0 和b = 1，而n=6 或步骤=6。但是，当它应该是= 4.288

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

    double f(double x)
    


        return exp(-x)/sqrt(1-x);
    

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

    double x0=f(a)+f(b);

    double x1=0,x2=0;
    double x=0;
    double h=(b-a)/(n);
    for(int i = 1 ; i <n;i++)

            x=a+(h*i);
            if(i%2==0)
            
                x2=x2+f(x);
            
            else
            
                x1=x1+f(x);
            


        
            return (h*(x0+2*x2+4*x1))/3;
;
    

     int main()
          cout<<endl;
               cout<<"The improper integral is: "<<simpson(0.00000009,1,6)<<" "<<endl;
               cout<<endl;
        

【问题讨论】：

您是否在 x=1 处进行评估？从h 和n 的定义看来你是。

它应该计算从下限 = 0 到上限 = 1 的积分，而 n=6

在 x=1 时，您的函数趋于无穷大。你需要保护它。

@macroland 但算法的目的不是计算不定积分吗？或其他可能的解决方案？

算法只是计算积分。只需给出一个公差，例如 1E-5 并从界限中减去。顺便说一句，n=6 似乎太少了。

【参考方案1】：

我认为您的问题是您正在评估的函数 f(double x) 未在 x = 1 处定义，因为它等于被零除的 exp(-1)/sqrt(1 - 1)，但您在第一行调用 f(b) simpson 当 b 为 1 时。因此，您返回无穷大。您可能希望将 b 设置为 0.9999。如果您使用更大的 n 值来更接近实际曲线（例如 n = 400），那么您将非常接近 @macroland 在 Wolfram Alpha 上找到的实际答案 1.076。

【讨论】：

我得到：不正确的积分是：2.92876

wolframalpha.com/input/…

谢谢@Evg。改变了答案。

谢谢，所以几乎不可能在 n=6 时得到准确的值？

这取决于函数@sanjep，但辛普森的方法在f(x)趋于无穷大的小范围内很可能不准确。如果你在 a = 0, b = 0.5 范围内对其进行测试，它应该接近完美。