插值型求积公式(C++)

插值型求积公式是一种用于计算定积分的数学方法. 它基于把未知的定积分分解为可求解的函数值的积分, 然后将函数拟合为多项式或其他函数, 最后使用拟合函数去近似定积分. 由于拟合的函数具有更好的变化性, 效果更好, 而改进的插值型求积公式又能更好地拟合变化的函数, 因此它们被认为是计算积分的有效方法.

说明

插值型求积公式要求输入解的定积分${\int }_a^{b}f(x)\mathrm{d}x$, 插值节点的个数$n$, 插值节点的选择方法, 以及插值函数的类型, 最终输出定积分的近似值. 这种求积方法导出的公式也被称为Newton-Cotes求积公式. 其算法可描述如下:

• 根据插值节点的选择方法, 在区间$[a,b]$上选择$n+1$个插值节点 $x_0,x_1,\cdots ,x_n$.
• 根据插值节点和插值函数的类型, 构造插值函数$P(x)$.
• 计算插值函数的积分${\int }_a^{b}P(x)\mathrm{d}x$.

我们可以通过插值公式计算出当$n$较小时的公式如下:
$${\int }_a^{b}f(x)\mathrm{d}x\approx \frac{b-a}{2}[f(a)+f(b)]$$
$${\int }_a^{b}f(x)\text{d}x\approx \frac{b-a}{6}[f(a)+4f(\frac{b+a}{2})+f(b)]$$
$${\int }_a^{b}f(x)\text{d}x\approx \frac{b-a}{8}\left[f(x_0)+3f(x_1)+3f(x_2)+f(x_3)\right]$$
$${\int }_a^{b}f(x)\text{d}x\approx \frac{b-a}{90}[7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)]$$

上述公式依次被称为梯形公式, Simpson公式, Newton公式和Cotes公式.

在实际应用中, 插值型求积公式通常会选择合适的插值节点和插值函数类型, 以使得近似值尽可能接近真实值. 此外, 还可以通过增加插值节点的个数$n$来提高近似精度. 但当插值节点增多时, 初始数据的误差引起计算结果的误差的增大, 即计算不稳定, 因此, 高阶Newton-Cotes公式不稳定. 所以当节点较多时, 需要寻求精度更高的求积公式.

代码实现

首先包含头文件:

#include <functional>
using std::function;

梯形公式:

/*
 * 梯形公式
 * f:被积函数
 * a:积分下限
 * b:积分上限
 */
double trapezoid(function<double(double)> f, const double &a, const double &b)
{
    if (a >= b)
        throw "区间上限小于区间下限！";
    double fa(f(a));
    if (isnan(fa))
        throw "积分下限处函数值无法计算！";
    double fb(f(b));
    if (isnan(fb))
        throw "积分上限处函数值无法计算！";
    return (b - a) * (fa + fb) / 2;
}

Simpson公式:

/*
 * Simpson公式
 * f:被积函数
 * a:积分下限
 * b:积分上限
 */
double Simpson(function<double(double)> f, const double &a, const double &b)
{
    if (a >= b)
        throw "区间上限小于区间下限！";
    double fa(f(a));
    if (isnan(fa))
        throw "积分下限处函数值无法计算！";
    double fb(f(b));
    if (isnan(fb))
        throw "积分上限处函数值无法计算！";
    double fm(f((a + b) / 2));
    if (isnan(fm))
        throw "积分区间中点处函数值无法计算！";
    return (b - a) * (fa + fb + 4 * fm) / 6;
}

Newton公式:

/*
 * Newton公式
 * f:被积函数
 * a:积分下限
 * b:积分上限
 */
double Newton(function<double(double)> f, const double &a, const double &b)
{
    if (a >= b)
        throw "区间上限小于区间下限！";
    double fa(f(a));
    if (isnan(fa))
        throw "积分下限处函数值无法计算！";
    double fb(f(b));
    if (isnan(fb))
        throw "积分上限处函数值无法计算！";
    double d0(b - a), d(d0 / 3), f1(f(a + d));
    if (isnan(f1))
        throw "积分区间左三等分点处函数值无法计算！";
    double f2(f(b - d));
    if (isnan(f2))
        throw "积分区间右三等分点处函数值无法计算！";
    return d0 * (fa + fb + (f1 + f2) * 3) / 8;
}

Cotes公式:

/*
 * Cotes公式
 * f:被积函数
 * a:积分下限
 * b:积分上限
 */
double Cotes(function<double(double)> f, const double &a, const double &b)
{
    if (a >= b)
        throw "区间上限小于区间下限！";
    double fa(f(a));
    if (isnan(fa))
        throw "积分下限处函数值无法计算！";
    double fb(f(b));
    if (isnan(fb))
        throw "积分上限处函数值无法计算！";
    double d0(b - a), d(d0 / 4), x(a + d), f1(f(x));
    if (isnan(f1))
        throw "积分区间左四等分点处函数值无法计算！";
    double fm(f(x += d));
    if (isnan(fm))
        throw "积分区间中点处函数值无法计算！";
    double f2(f(x += d));
    if (isnan(f2))
        throw "积分区间右四等分点处函数值无法计算！";
    return d0 * (7 * (fa + fb) + 32 * (f1 + f2) + 12 * fm) / 90;
}