梯形规则的数值积分c++

Numerical Integration with Trapezoidal rule c++

本文关键字:数值积分 c++ 规则      更新时间:2023-10-16

我正在尝试实现一个梯形规则,该规则利用以前的函数求值来避免冗余计算。有两件事:a)计算结果不收敛,我有点不确定为什么。我将公布为什么我认为如果需要的话,算法应该产生收敛性,以及b)do while循环终止于n=8,我也无法计算出这一点;它应该一直运行到n>128?(n是子区间的数量)我的代码如下。提前感谢!

void NestedTrap(int n) //Trapezoidal with reuse of function evaluations
{
    double a,b; //interval end points
    double x[n+1]; //equally spaced nodes
    double c[n]; //midpoints
    double T; //Initial integral evaluation
    double T2; //Evaluation with reuse of previous function evaluations
    double h, h2; //step sizes for T and T2
    double temp1, temp2;
    std::cout <<"Enter interval end points (lesser first; enter 999 for pi & 999.2 for pi/2 & 999.4 for pi/4): ";
    std::cin >> a >> b;
    if (b == 999)
    {
        b = M_PI;
    }
    if (a == 999)
    {
        a = M_PI;
    }
    if (b == 999.4)
    {
        b = M_PI/4;
    }
    if (a == 999.4)
    {
        a = M_PI/4;
    }
    if (b == 999.2)
    {
        b = M_PI/2;
    }
    if (a == 999.2)
    {
        a = M_PI/2;
    }
    h = (b-a)/n;
    T = 0;
    temp1 = 0;
    temp2 = 0;
    for (int i=0; i<=n; i++)
    {
        x[i] = 0;
    }

    for (int i=0; i<n; i++)
    {
        x[i+1] = x[i] + h;
    }
    for (int i=1; i<n; i++)
    {
        temp1 += I1(x[i]);
    }
    T = (h/2)*exp(x[0]) + (h/2)*exp(x[n]) + (h*temp1);
    std::cout << "T_" << n <<": " << T << std::endl;
    do
    {
    temp2 = 0;
    n = 2*n;
    h2 = (b-a)/(n);
    for (int i=0; i<n; i++)
    {
        c[i] = 0;
    }
    for (int i=1; i<=n; i++)
    {
        c[i] = a + h2*(i-0.5);
        //std::cout << c[i] << std::endl;
    }
    for (int i=0; i<n; i++)
    {
        temp2 += exp(c[i]);
    }
    T2 = (T/2) + h2*temp2;
    std::cout << "T_" << n <<": " << T2 << std::endl;
    T = T2;
    }   while (n <= 128);
}

您在此处创建大小为n的数组

double x[n+1]; //equally spaced nodes
double c[n]; //midpoints

(注意这不是有效的c++)

然后在此处增加n

n = 2*n;

然后你在这里写过数组的末尾:

for (int i=0; i<n; i++)
{
    c[i] = 0;
}

这会导致未定义的行为(可能会覆盖其他一些变量)