龙贝格算法C语言实现(数值分析经典算法)

算法理论

注意点

需要明确区别变量k和m：
k 表示区间[a,b]的二分数，等价于区间[a,b]被分为2^k 个小区间，等价于每个小区间的步长h = (a-b)/2^k。
m 表示理查森外推加速方法的加速次数，在此处的梯形值可直接由查理森外推算法的公式进行行递归而得。但需要注意的是，在进行行递归的时候一定要有初始值（即在二分数为k 时的梯形初始值）。

算法应用实例

实例

定义函数接口：

double Integral(double a, double b, double(*f)(double x, double y, double z), double TOL, double l, double t)

在接口定义中：a、b 分别为定积分的上、下界，f 是积分核函数，其中x是积分哑元，y、z 是本题目定义的特殊参数，分别对应$y(x)=l\ast sin(tx)$中的 $l$ 和 $t$ 的值。另外需要注意，在本实例情景下，$y(x)=l\ast sin(tx)$的输出单位为厘米，而最终的输出的长度值要求以米为单位，需要进行单位换算。

解读实例

实例即可转化为数学问题：已知积分${\int }_a^{b}f(x)dx$,用龙贝格求积算法，并结合梯形公式求积分在精度为$TOL$时的近似值。其中，需要注意的是，$f(x)=y(x)=l\ast sin(tx)$，且 $l$ 与 $t$ 为给定的参数。

代码思路分析

整个代码实现分为四个部分，分别为
1、弧长积分函数

double f0( double x, double l, double t )

2、梯形公式函数

double T(double a, double b, double h, double (*f)(double x, double y, double z), double l, double t)

3、龙贝格算法实现函数

double Integral(double a, double b, double (*f)(double x, double y, double z), double TOL, double l, double t)

4、主函数

int main()

其中核心部分为梯形公式函数 $T$ 和龙贝格算法的实现函数 $Integral$ 。

(1)梯形公式函数 T 部分：

此部分只需根据复合梯形公式
$T=\frac{h}{2}[f(a)+2\ast {\sum }_{k=1}^{n-1}f(x_k)+f(b)]$
并并调用弧长积分函数 $f0$ 即可编写出，其中需要注意的是，由于 $n$ 与 $h$ 之间有关系 $h=\frac{(b-a)}{n}$，故可以在T 函数的参数传入部分只需传入 $n$ 与 $h$ 中的任意一个与其他参数组合即可，又由于后续龙贝格算法实现函数Integral 中通过 $h$ 作为参数来调用梯形公式函数 $T$ 更方便，所以在该部分采用 $h$ 作为参数之一。

(2)龙贝格算法实现 Integral 部分

此部分需要调用弧长积分函数 $f0$、梯形公式函数 $T$ 以及采用递归的方法实现。分为以下步骤：

Step1:

计算二分数 $k=0$ 时的梯形值 $T0$，这一步通过调用T 函数即可实现。

Step2:

计算二分数 $k=1$ 时的初始梯形值 $T1\_0$（通过调用 $T$ 函数实现），并计算经过理查森外推加速方法第一次加速后的梯形值 $T1$（通过调用 $T$ 函数以及理查森外推加速方法公式实现）。

Step3:

建立循环，当 $|T1-T0|<TOL$ 不成立时便会一直处于循环当中。在循环里通过进行 $T0=T1$ 以及重置
$T1\_0$ 的操作来实现变量的重利用和递归，并通过建立for 循环实现k 次外推，最终返回满足 $|T1-T0|<TOL$ 的积分近似值（返回时应进行单位换算）。

代码

#include<stdio.h>
#include<math.h>

double f0( double x, double l, double t )
{ /* 弧长积分函数 */
    return sqrt(1.0+l*l*t*t*cos(t*x)*cos(t*x));
}

double T(double a, double b, double h, double (*f)(double x, double y, double z), double l, double t)
{ /*梯形公式函数*/
    double x; //分点
    double sum = 0; //求除a、b外的其他分点对应的函数总和
    int n = (b-a)/h;
    int k;
    for(k=1; k<n; k++)
    {
        x = a+k*h;
        sum += f(x,l,t);
    }
    double I = h/2*(f(a,l,t)+2*sum+f(b,l,t));
    return I;
}

double Integral(double a, double b, double (*f)(double x, double y, double z), double TOL, double l, double t)
{ /*龙贝格算法实现函数*/
    int k = 0;
    double h;
    h = (b-a)/pow(2,k); //h为小区间长度
    double T0 = T(a,b,h,f0,l,t);
    
    k = 1;
    h = (b-a)/pow(2,k);
    double T1_0 = T(a,b,h,f0,l,t);
    double T1 = (pow(4,k)/(pow(4,k)-1)*T(a,b,h/2,f0,l,t)) - (1/(pow(4,k))*T1_0);

    int i;
    while(fabs(T1-T0)>=TOL)
    {
        k++;
        h = (b-a)/pow(2,k);
        T0 = T1;
        T1_0 = T(a,b,h,f0,l,t);
        for(i=1; i<=k; i++)
        {
            T1 = (pow(4,i)/(pow(4,i)-1)*T(a,b,h/2,f0,l,t)) - (1/(pow(4,i))*T1_0);
            T1_0 = T1;
        }
    }
    h = (b-a)/pow(2,k);
    return T(a,b,h,f0,l,t)/100;
}

int main()
{
    double a=0.0, b, TOL=0.005, l, t;
    while (scanf("%lf %lf %lf", &l, &b, &t) != EOF)
        printf("%.2f\n", Integral(a, b, f0, TOL, l, t));
    return 0;
}

测试

Input：2 100 1
Output：1.68 

初来乍到，如有不足，欢迎大家批评指正~