Рекурсия для функции аппроксимации интегрирования

#c #math

#c #математика

Вопрос:

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

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

 //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 , поэтому все улучшения приветствуются.

Комментарии:

1. Где определен currentLevel? Я не вижу этого нигде в вашем коде. Если это глобальная переменная, вы делаете что-то не так.

2. Да, у меня есть currentLevel, определенный как глобальная переменная. Спасибо за ваш ответ, приведенный ниже. Я рассмотрю это более тщательно.

Ответ №1:

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

Комментарии:

1. Я понимаю. Глупая математическая ошибка при вычислении интервалов. Почему плохо увеличивать глобальную переменную? Кажется, это приведет к тому же эффекту. Большое вам спасибо за вашу помощь.

2. Я вижу, что помогло использование вспомогательной функции, так что мне не нужно, чтобы currentLevel был глобальной переменной. Теперь он сошелся на правильном значении, однако я не понимаю, почему это помогло.

3. Теперь я понимаю, почему я не хочу, чтобы currentLevel был глобальной переменной. Наличие его в качестве глобальной переменной полностью противоречит концептуальной цели адаптивного алгоритма. Мне нужно, чтобы это зависело от текущего уровня рекурсии функции. Спасибо.

4. @kienjakenobi: Да, именно по этой причине. В более сложных сценариях, если вы оставите это как глобальную переменную, станет еще более очевидным, почему эта схема сломается. Желаю удачи в вашем будущем программировании на C !