C语言 递归查找二阶导数

Question

>我正在尝试用 C 编写一个程序，该程序可以找到给出 x 的 f（x） 二阶导数的数值。当我使用我的函数找到一阶导数时，一切似乎都正常，但递归二阶导数部分总是给我 0.0。这是我的函数：

double derivative(double (*f)(double), double x0, int order) {
    if(order == 1){
        const double delta = 1.0e-6;
        double x1 = x0-delta;
        double x2 = x0+delta;
        double y1 = f(x1);
        double y2 = f(x2);
        return (y2 - y1) / (x2 - x1);
    } else if(order == 2) {
        const double delta = 1.0e-6;
        double x1 = derivative(f, x0, 1)-delta;
        double x2 = derivative(f, x0, 1)+delta;
        double y1 = derivative(f, f(x1), 1)-delta;
        double y2 = derivative(f, f(x2), 1)+delta;
        return (y2 - y1) / (x2 - x1);
    } else {
        printf("order too high error n");
        return 0;
    }
}

我

理解为什么它不能按预期工作——在这里，我试图通过基于函数 f 和 f'（x） 的值来找到 f 的二阶导数。但是我已经为此工作了几个小时，我无法想出合适的算法来做到这一点。

任何帮助将不胜感激。提前感谢！

Answer 1

我将您的代码更改为以下代码，它按预期工作。

double derivative(double (*f)(double), double x0, int order)
{
         const double delta = 1.0e-6;
         if (order == 1) {
                  double x1 = x0 - delta;
                  double x2 = x0 + delta;
                  double y1 = f(x1);
                  double y2 = f(x2);
                  return (y2 - y1) / (x2 - x1);
          } else if(order == 2) {
                  double x1 = x0 - delta;
                  double x2 = x0 + delta;
                  double y1 = derivative(f, x1, 1);
                  double y2 = derivative(f, x2, 1);
                  return (y2 - y1) / (x2 - x1);
          } else {
                  printf("order too high error n");
                  return 0;
          }
  }

Answer 2

我会像这样实现它：

double derivative(double (*f)(double), double x0, int order)
{
         const double delta = 1.0e-6;
         double x1 = x0 - delta;
         double x2 = x0 + delta;
         if (order == 1) {
                  double y1 = f(x1);
                  double y2 = f(x2);
                  return (y2 - y1) / (x2 - x1);
          } else {
                  double y1 = derivative(f, x1, order - 1);
                  double y2 = derivative(f, x2, order - 1);
                  return (y2 - y1) / (x2 - x1);
          }
 }

因为这样就不需要"订单太高的错误"。

Answer 3

Parham Alvani 已经回答了你的问题，并指出了你的函数中的逻辑错误。所以这不是一个真正的答案，只是一些额外的评论。

您可以使函数完全递归，而无需单独处理除基本情况之外的任何情况：

double derivative(double (*f)(double), double x, int order)
{
    const double eps = 1.0e-3;
    if (order == 0) return f(x);
    double y1 = derivative(f, x - eps, order - 1);
    double y2 = derivative(f, x + eps, order - 1);
    return (y2 - y1) / (2 * eps);
}

但是，对于高阶导数，您最终会多次计算相同的函数值，因此您可以递归地计算导数并对关系进行硬编码：

  Δf(x) / Δx = (f(x + ε) - f(x - ε)) / (2·ε)
Δ²f(x) / Δx² = (f(x + 2·ε) - 2·f(x) + f(x - 2·ε)) / (4·ε²)
Δ³f(x) / Δx³ = (f(x + 3·ε) - 3·f(x + ε) + 3·f(x - ε) - f(x - 3·ε)) / (8·ε³)

等等。然而，对于高阶导数，选择ε值变得棘手。当计算的函数值之间的差异很小时，选择太小将导致错误。