積分近似関数の再帰

Question

適応台形法を用いて近似を近似しようとしています。

Iが粗い積分近似有する：元気積分近似を有する

//Approximates the integral of f across the interval [a,b] 
double coarse_app(double(*f)(double x), double a, double b) { 
    return (b - a) * (f(a) + f(b))/2.0; 
} 

を：

//Approximates the integral of f across the interval [a,b] 
double fine_app(double(*f)(double x), double a, double b) { 
    double m = (a + b)/2.0; 
    return (b - a)/4.0 * (f(a) + 2.0 * f(m) + f(b)); 
} 

これはれるまで所定の間隔の減少する部分を横切って近似値を合計することによって適応作られますか再帰レベルが高すぎるか、または粗いおよび細かい近似がお互いに非常に接近している：

//Adaptively approximates the integral of f across the interval [a,b] with 
// tolerance tol. 
double trap(double(*f)(double x), double a, double b, double tol) { 
    double q = fine_app(f, a, b); 
    double r = coarse_app(f, a, b); 
    if ((currentLevel >= minLevel) && (abs(q - r) <= 3.0 * tol)) { 
     return q; 
    } else if (currentLevel >= maxLevel) { 
     return q; 
    } else { 
     ++currentLevel; 
     return (trap(f, a, b/2.0, tol/2.0) + trap(f, a + (b/2.0), b, tol/2.0)); 
    } 
} 

積分を手動でセクションに分割し、それにfine_appを使用して計算すると、非常に良い近似が得られます。しかし、私のためにこれを行う必要があるトラップ関数を使用すると、私の結果はすべてあまりにも小さいです。

たとえば、trap（square、0、2.0、1.0e-2）は出力0.0424107を返します。ここで、二乗関数はx^2として定義されています。しかし、出力は約2.667になるはずです。これは、全体の間隔でfine_appを1回実行するよりもはるかに悪く、値は3です。

概念的には、正しく実装されていると思いますが、C++の再帰については何かがありません。それを期待する。

初めてC++でプログラミングするので、すべての改善が歓迎されます。

Answer 1

私はcurrentLevelが別の場所で定義されていると仮定しています。あなたはそれをしたくありません。また、中点を間違って計算します。

取る= 3、B = 5：

[a, b/2.0] = [3, 2.5] 
[a + b/2.0, b] = 2.5, 3] 

正しい点

がなければならない[3、4]及び[4、5]

コードは次のようになります

double trap(double(*f)(double x), double a, double b, double tol, int currentLevel) { 
    double q = fine_app(f, a, b); 
    double r = coarse_app(f, a, b); 
    if ((currentLevel >= minLevel) && (abs(q - r) <= 3.0 * tol)) { 
     return q; 
    } else if (currentLevel >= maxLevel) { 
     return q; 
    } else { 
     ++currentLevel; 
     return (trap(f, a, (a + b)/2.0, tol/2, currentLevel) + trap(f, (a + b)/2.0, b, tol/2, currentLevel)); 
    } 
} 

あなたがcurrentLevelを指定する必要はありませんので、あなたはヘルパー関数を追加することができます。

double integrate(double (*f)(double x), double a, double b, double tol) 
{ 
    return trap(f, a, b, tol, 1); 
} 

私がこれをintegrate(square, 0, 2, 0.01)と呼んだ場合、答えは2.6875になります。つまり、正確な結果である8/3 = 2.6666...7に収束するにはさらに低い許容差が必要です。 Simpson's methodのエラー条件を使用して、正確なエラー境界を確認できます。