递归求二阶导数
Recursively find the second derivative
我正在尝试用 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 的二阶导数。但是我已经为此工作了几个小时,但我想不出一个合适的算法来做到这一点。
如有任何帮助,我们将不胜感激。提前致谢!
我将您的代码更改为以下代码,它按预期工作。
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;
}
}
我会这样实现:
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);
}
}
因为那样就不需要 "order too high error".
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·ε³)
等等。然而,对于高阶导数,选择 ε
值会变得棘手。选得太小会在求值函数值相差较小时产生误差
我正在尝试用 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 的二阶导数。但是我已经为此工作了几个小时,但我想不出一个合适的算法来做到这一点。
如有任何帮助,我们将不胜感激。提前致谢!
我将您的代码更改为以下代码,它按预期工作。
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;
}
}
我会这样实现:
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);
}
}
因为那样就不需要 "order too high error".
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·ε³)
等等。然而,对于高阶导数,选择 ε
值会变得棘手。选得太小会在求值函数值相差较小时产生误差