我實作了一個 C 代碼來數值求解點 x_0 中函式的 n 次導數:
double n_derivative( double ( *f )( double ), double x_0, int n )
{
if( n == 0 ) return f( x_0 );
else
{
const double h = pow( __DBL_EPSILON__, 1/3 );
double x_1 = x_0 - h;
double x_2 = x_0 h;
double first_term = n_derivative( f, x_2, n - 1 );
double second_term = n_derivative( f, x_1, n - 1);
return ( first_term - second_term ) / ( 2*h );
}
}
我想知道這是否適合您的一個很好的實作,或者是否有辦法更好地用 C 撰寫它。問題是我注意到 n 次導數對于 n 大于 3 的值會發散。你知道如何解決這個問題嗎?
uj5u.com熱心網友回復:
這不是一個好的實作
至少這些問題。
整數數學
使用 FP 數學作為1/3零。
1/3 --> 1.0/3
使用立方根優化 n==1
但不一定是其他n的。 @尤金
錯誤的ε
下面的代碼僅對|x_0|1.0 左右有用。x_0大時,x_0 - h可等于x_0。x_0小時,可x_0 - h相等。-h
OP 的 /- some epsilon 適用于定點,但double它是浮點數。
// Bad
const double h = pow( __DBL_EPSILON__, 1.0/3 );
double x_1 = x_0 - h;
需要相對縮放。
#define EPS cbrt(DBL_EPSILON) // TBD code to well select this
if (fabs(x_0) >= DBL_MIN && isfinite(x_0)) {
double x_1 = x_0*(1.0 - EP3);
double x_2 = x_0*(1.0 EPS);
double h2 = x_2 - x_1;
...
} else {
TBD_Code for special cases
}
無效的代碼
f是double ( *f )( int, double ),但呼叫是f( x_0 )
次要:混淆名稱
為什么first_termwithx_2和second_termwith x_1?
轉載請註明出處,本文鏈接:https://www.uj5u.com/gongcheng/417537.html
標籤:
