0

所以我有一个程序在函数 x^2 + y^2 < 1 上实现自适应 2D 梯形规则，但似乎递归不起作用——这里的程序是（工作）1D 的修改形式梯形方法，所以我不确定代码在哪里发生故障，它应该返回 PI：

double trapezoidal(d_fp_d f,
                   double a, double b,
                   double c, double d) { //helper function
    return 0.25*(b-a)*(d-c)*
    (f(a, c)+f(a, d) +
     f(b, c)+f(b, d));
}

double atrap( double a, double b, double c, double d, d_fp_d f, double tol )
 {// helper function

 return atrap1(a, b, c, d, f, tol );
 }
double atrap1( double a, double b, double c, double d, d_fp_d f, double tol)
{
 //implements 2D trap rule 
    static int level = 0;
    const static int minLevel = 4;
    const static int maxLevel = 30;
    ++level;
    double m1 = (a + b)/2.0;
    double m2 = (c + d)/2.0;
    double coarse = trapezoidal(f,a,b,c,d);
    double fine = 
      trapezoidal(f, a, m1, c, m2)
    + trapezoidal(f, a, m1, m2, d)
    + trapezoidal(f, m1, b, c, m2)
    + trapezoidal(f, m1, b, m2, d);
    ++fnEvals;
    if( level< minLevel
       || ( abs( fine - coarse ) > 3.0*tol && level < maxLevel ) ){

            fine =  atrap1( a, m1, c, m2, f,tol/4.0)
            + atrap1( a, m1, m2, d, f, tol/4.0)
            + atrap1(m1, b, c, m2, f, tol/4.0)
            + atrap1(m1, b, m2, d, f,tol/4.0);

        }

    --level;
    return fine;
}

函数由下式给出

double ucircle( double x, double y)
{
    return x*x + y*y < 1 ? 1.0 : 0.0;
}

我的主要功能是

int main()
{
   double a, b, c, d;
    cout << "Enter a: " <<endl;
    cin >> a;
    cout << "Enter b: " <<endl;
    cin >> b;
    cout << "Enter c: " <<endl;
    cin >> c;
    cout << "Enter d: " <<endl;
    cin >> d;

    cout << "The approximate integral is: " << atrap( a, b, c, d, ucircle, 1.0e-5) << endl;

    return 0;
}
4

2 回答 2

1

它实际上不会永远运行，但它实际上运行了很长时间，您认为它会永远运行，这就是原因：在第一次运行时level，函数输入你的if，它调用自己 4 次，现在考虑第一次：它也进入if并再调用自身 4 次并继续......对于正确选择的输入（如您指定的输入），条件abs(fine - coarse)始终true是唯一可以阻止流进入的东西if是level会增加然后减少所以你的函数几乎会被调用4^30，这真的是一个很大的数字，你在一两个小时内看不到它的结束！

于 2012-10-21T01:40:57.807 回答
0

就像 BigBoss 已经写的那样，你的程序应该完成，它只需要很长时间，因为 30 次递归意味着4^30对 的函数调用atrap1，即1152921504606846976. 让这个数字沉入其中。

这里还有一些需要考虑的事情：

  • 您可能想在“中断条件”中使用fabs而不是。abs（我认为您应该收到有关整数转换的警告 - 或类似的东西 - 为此）abs可能会返回不可预测的值float或double参数。非常高的价值。

  • tol似乎是一个代表目标精度值的变量。但是，随着每次递归，您会进一步提高此目标精度。在第 10 次递归时，它已经是1E-11. 不确定这是有意的。随便什么tol意思。

    您可能不希望递归调用中的/4.0（.0顺便说一句是多余的）。

  • 您确实通过优化来编译它，对吗？

  • trapezoidal, minLevel,maxLevel可以是宏。

  • level由于是静态的，您的函数不喜欢线程执行。您应该将其作为atrap1.

于 2012-10-22T09:11:19.607 回答