数值积分中,插值求积公式可以近似求积函数。对于插值函数,可以用梯形公式来近似,为了满足一定的精度,需要将求积区间分成很多个小的区间,分别用于梯形公式,但是由于某些时候步长不能很好的控制,所以需要变步长。本算法给出变步长求积公式,该算法的具体形式可以在数值分析的书籍中找到,下面直接给出算法的c++代码;
/*变步长梯形求积法
*该函数double tarp_form(double a,double b,double eps,double (*f)(double)
*double a:积分区间左端点
*double b:几根区间右端点
*double eps:误差精度
*double* f(double):积分函数
*函数值返回double类型
*/
#include<iostream>
#include<cmath>
double trap_form(double a, double b, double eps, double (*f)(double))
{
int n, k;
double fa, fb, h, t1, p, s, x, t;
fa = (*f)(a); //计算左端点的值
fb = (*f)(b); //计算右端点的值
n = 1;
h = b - a; //初始步长
t1 = h*(fa + fb) / 2.0; //第一次积分
p = eps + 1.0;
while (p >= eps)
{
s = 0.0;
for (k = 0; k <= n - 1; k++)
{
x = a + (k + 0.5)*h;
s = s + (*f)(x);
}
t = (t1 + h*s) / 2.0;
p = fabs(t1 - t); //前后积分之差
//继续增加结点
t1 = t;
n = 2 * n;
h = h / 2.0;
}
return t;
}
//积分函数y=exp(-x^2)
double func1(double x)
{
return exp(-x*x);
}
//积分函数y=1/(1+x^2)
double func2(double x)
{
return 1 / (1 + x*x);
}
//func3=sinx/x
double func3(double x)
{
return sin(x) / x;
}
int main()
{
double a = 0.0;
double b = 1.0000000;
double eps = 0.000001;
double t = trap_form(a, b, eps, func2);
std::cout << "t= " << t << std::endl;
t = trap_form(a, b, eps, func1);
std::cout << "func1's t= " << t << std::endl;
//由于func3的分母是x,a计算的时候不能为0,所以我们只能无限接近于0,否则会计算不出值
a = 0.000000001;
t = trap_form(a, b, eps, func3);
std::cout << "func3's value " << t << std::endl;
return 0;
}
下面是运算结果: