复合辛普森规则无限输出
Composite simpson rule infinity output
算法:
我有以下代码使用复合辛普森规则计算不正确的积分,我正在尝试计算积分 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;
}
我认为你的问题是你正在评估的函数 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。
算法:
我有以下代码使用复合辛普森规则计算不正确的积分,我正在尝试计算积分 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;
}
我认为你的问题是你正在评估的函数 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。