Блог

Я пытаюсь аппроксимировать интегралы, используя адаптивное трапециевидное правило.

У меня есть грубое интегральное приближение:

 //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) amp;amp; (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, которая должна делать это за меня, все мои результаты слишком малы.

Например, trap(square, 0, 2.0, 1.0e-2) выдает результат 0.0424107, где квадратная функция определяется как x ^ 2. Однако результат должен быть около 2.667. Это намного хуже, чем выполнение одного запуска fine_app на всем интервале, что дает значение 3.

Концептуально, я считаю, что я реализовал ее правильно, но в рекурсии C есть что-то, что не делает то, что я ожидаю.

Впервые программирую на C , поэтому все улучшения приветствуются.

Я предполагаю, что у вас есть currentLevel, определенный где-то еще. Вы не хотите этого делать. Вы также неправильно вычисляете свои средние точки.

Возьмем a = 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) amp;amp; (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 . Вы можете проверить точную ошибку, связанную с этим, используя условия ошибки для метода Симпсона.