（计算机数值分析）龙贝格算法C语言版

#1024程序员节｜征文#

计算机数值分析课上作业。

在讲完了复化梯形公式T，复化辛甫生公式S和科特斯公式C之后，观察他们的公式发现存在一个关系。即：

(I-T2n)/I-Tn = 1/4,(I-S2n)/(I-Sn) = 1/16,(I-C2n)/(I-Cn) = 1/64

于是得到以下关系：Sn = 3/4T2n - 1/3Tn
                                   Cn = 16/15S2n - 1/15Sn

由此得到龙贝格公式：Rn = 64/63C2n - 1/63Cn

运用龙贝格公式可以将收敛速度缓慢的梯形公式迅速收敛。

一下就例题展示。

如果仅使用复化梯形公式计算积分值，需要用变步长方法二分10次才能得到结果。但是如果我们使用龙贝格公式的话，能加速计算，大大降低时间和运算成本。

然后我们仅需要使用二分3次的数据，通过3次加速就获得了二分10次才能求得的结果。

以下是c语言源代码呈现：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
long double Value(double x)
{
    return sin(x)/x;
}

int main(void) {
    int a=1, b=1,k=1;
    double h=1,x=1;
    long double e=1,S=1,S1=1,S2=1,T1=1,T2=1,C1=1,C2=1,R1=1,R2=1;
    printf("请输入a,b,e:");
    scanf("%d%d%Lf", &a, &b, &e);
    h = b - a;
    T1 = (1 + Value(b)) * (h / 2);
    printf("T1的值为：%.7Lf\n",T1);
    k = 1;
    S = 0;
    x = a + h / 2;
    while (x < b)
    {
    S += Value(x);
    x += h;
    }
    T2 = T1/2+(h/2)*S;
    S2 = T2 + 1/3*(T2-T1);
    if(k==1){
        k += 1;
        h = h/2;
        T1 = T2;
        S1 = S2;
        S = 0;
        x = a + h / 2;
        while (x < b)
        {
            S += Value(x);
            x += h;
        }
        T2 = T1/2+(h/2)*S;
        S2 = T2 + (T2-T1)/3;
    }
    C2 = S2 + (S2 - S1)/15;
    if(k==2)
    {
        C1 = C2;
        k += 1;
        h = h/2;
        T1 = T2;
        S1 = S2;
        S = 0;
        x = a + h / 2;
        while (x < b)
        {
            S += Value(x);
            x += h;
        }
        T2 = T1/2+(h/2)*S;
        S2 = T2 + (T2-T1)/3;
        C2 = S2 + (S2-S1)/15;
    }
    R2 = C2 + (C2 - C1)/63;
    if(k==3)
    {
        R1 = R2;
        C1 = C2;
        k += 1;
        h = h/2;
        T1 = T2;
        S1 = S2;
        S = 0;
        x = a + h / 2;
        while (x < b)
        {
            S += Value(x);
            x += h;
        }
        T2 = T1/2+(h/2)*S;
        S2 = T2 + (T2-T1)/3;
        C2 = S2 + (S2-S1)/15;
        R2 = C2 + (C2 - C1)/63;
    }
    while(abs(R2 - R1) >= e)
    {
        R1 = R2;
        C1 = C2;
        k += 1;
        h = h/2;
        T1 = T2;
        S1 = S2;
        S = 0;
        x = a + h / 2;
        while (x < b)
        {
            S += Value(x);
            x += h;
        }
        T2 = T1/2+(h/2)*S;
        S2 = T2 + (T2-T1)/3;
        C2 = S2 + (S2-S1)/15;
        R2 = C2 + (C2 - C1)/63;
    }
    printf("R2的值为：%.7Lf",R2);
    return 0;
}

运行之后输入：0 1 1
出现以下结果：

对照给出的数据表：
​​​​​​​

发现计算出的值恰好是数值表中给出的值。