龙贝格(Romberg)算法的应用

512 阅读3分钟

实验要求:

下面图1中的塑料雨蓬材料是由图2中所示的长方形平板塑料材料压制而成。

image.png 图一 image.png 图二

已知图1的横截面曲线形状满足函数,则给定了雨蓬的长度后,要求需要平板原材料的长度。

函数接口定义:

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是本题目定义的特殊参数,分别对应中的 lt;TOL是要求积分达到的精度;l和t传入裁判输入的参数 lt 的值。

image.png 另注意:

image.png的单位是厘米,输出的长度值要求以米为单位。

实验代码

#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 Function_add(double a ,double b ,double (*f)(double x, double y, double z),double h,double l,double t)

{

    double sum=0;

    double i=a+h;

    while(i<b)

    {

       sum=sum+f(i,l,t);

       i=i+h;

    }

    return sum;

}

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

{

    double h;

    double Tk[100][100];

    double flog,ans;

    int tg=1;

    h=b-a;

    Tk[0][0]=(h/2)*(f(a,l,t)+f(b,l,t));

    int i,j;

    for(i=1;i<100;i++)

    {

        h=h/2;

        Tk[i][0]=0.5*Tk[i-1][0]+0.5*h*Function_add(a,b,f,h,l,t);/*用于求算T[k][0]*/

        for(j=1;j<=i;j++)

        { /*按行进行外推*/

            Tk[i][j]=(pow(4,j)/(pow(4,j)-1))*Tk[i][j-1]-(1/(pow(4,j)-1))*Tk[i-1][j-1];

            flog=Tk[i][j]-Tk[i-1][j-1];

            ans=Tk[i][j]/100;

            if(fabs(flog)<=TOL)

            {

                tg=0;

                break;

            }

        }

        if(tg==0)

        {

            break;

        }

    }

    return ans;

}

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;

}

程序框图:

image.pngimage.png

实验心得和感想:

本次实验主要的目的是通过复合梯形公式结合理查德外推加速算法来求函数的积分。复合梯形公式的主要思想是将积分区间不断对分,如何每个区间的积分面积通过这个区间的梯形面积来近似。随着积分区间的不断细分,梯形的面积也不断的接近区间积分的值,当前一个算出的值与后一个算出的值的差小于给出的误差范围时,此时复合梯形公式算出的值可以近似的等于目标函数在积分区间的积分值。而理查德外推算法主要是通过:

image.png

image.png

image.png 的规律来实现对复合梯形公式的加速。在编程过程中应当注意的问题:其实本次实验代码的编写难度相对于前两次实验来说并不大,因为只要按着老师上课给的算法编写就可以了。但是在编程过程中特变应该注意下标的问题,应该稍微有一个下标出错就有可能导致结果出错并且这种错误很难检查出来,所以在实验过程中如果出现这种错误可以在关键节点加入打印语句来帮助检查。其次在编程过程中要特别注意h的取值看是带h还是代h/2。我在编程过程中就犯了这样的错误导致我的结果一值是3点多和标准结果相差很大,最后才发现是h代成了h/2.