数值微分

如何计算在无穷远处涉及指数和奇点的函数的数值二阶导数。 不幸的是,Ridder的方法在“C中的数值配方”中提供的数值导数只能计算一阶导数(它需要预先对函数进行解析表达式。)此外,我尝试了Chebyshev近似并在之后区分函数,但给出的值是方式关闭实际值。 我也尝试过在数学论文中提供的一些有限差分算法,但它们也容易出错。 函数是e ^(x / 2)/ x ^ 2。 对此事我感激不尽。

提前致谢

最新编辑:问题解决了C ++中提供的FADBAD库做得非常好。 它们可以通过http://www.fadbad.com/fadbad.html获得

编辑:

// The compilation command used is given below // gcc Q3.c nrutil.c DFRIDR.c -lm -o Q3 #include  #include  #include "nr.h" #define LIM1 20.0 #define a -5.0 #define b 5.0 #define pre 100.0 // This defines the pre /* This file calculates the func at given points, makes a * plot. It also calculates the maximum and minimum of the func * at given points and its first and second numerical derivative. */ float func(float x) { return exp(x / 2) / pow(x, 2); } int main(void) { FILE *fp = fopen("Q3data.dat", "w+"), *fp2 = fopen("Q3results.dat", "w+"); int i; // Declaring our loop variable float x, y, min, max, err, nd1, nd2; // Define the initial value of the func to be the minimum min = func(0); for(i = 0; x < LIM1 ; i++) { x = i / pre; // There is a singularity at x = 0 y = func(x); if(y < min) min = y; fprintf(fp, "%f \t %f \n", x, y); } fprintf(fp, "\n\n"); max = 0; for(i = 0, x = a; x  max) max = y; } fprintf(fp2, "The minimum value of f(x) is %f when x is between 0 and 20. \n", min); fprintf(fp2, "The maximum value of f(x) is %f when x is between -5 and 5. \n", max); fclose(fp); fclose(fp2); return 0; } 

编辑:切比雪夫

 // The compilation command used is given below //gcc Q3.c nrutil.c CHEBEV.c CHEBFT.c CHDER.c -lm -o Q3 #include  #include  #include "nr.h" #define NVAL 150 // Degree of Chebyshev polynomial #define LIM1 20.0 #define a -5.0 #define b 5.0 #define pre 100.0 // This defines the pre /* This file calculates the func at given points, makes a * plot. It also calculates the maximum and minimum of the func * at given points and its first and second numerical derivative. */ float func(float x) { return exp(x / 2) / pow(x, 2); } int main(void) { FILE *fp = fopen("Q3data.dat", "w+"), *fp2 = fopen("Q3results.dat", "w+"); int i; // Declaring our loop variable float x, y, min, max; float nd1, nd2, c[NVAL], cder[NVAL], cder2[NVAL]; // Define the initial value of the func to be the minimum min = func(0); for(i = 0; x < LIM1 ; i++) { x = i / pre; // There is a singularity at x = 0 y = func(x); if(y < min) min = y; fprintf(fp, "%f \t %f \n", x, y); } fprintf(fp, "\n\n"); max = 0; // We make a Chebyshev approximation to our function our interval of interest // The purpose is to calculate the derivatives easily chebft(a,b,c,NVAL,func); //Evaluate the derivatives chder(a,b,c,cder,NVAL); // First order derivative chder(a,b,cder,cder2,NVAL); // Second order derivative for(i = 0, x = a; x  max) max = y; } fprintf(fp2, "The minimum value of f(x) is %f when x is between 0 and 20. \n", min); fprintf(fp2, "The maximum value of f(x) is %f when x is between -5 and 5. \n", max); fclose(fp); fclose(fp2); return 0; } 

该function是可区分的,因此使用数值方法可能不是最好的。 二阶导数是:

6 * exp(x / 2)/(x ^ 4)-2 * exp(x / 2)/ x ^ 3 + exp(x / 2)/(4 * x ^ 2)

当然可以简化上述内容以加速计算。 编辑:第一次有原始配方错误。

如果您想要100%数值方法,那么请查看公共样条插值的数字配方 (Charter 3.3)。 它会在任何位置为您提供第二个衍生物。

使用xy值调用spline()以返回y2的第二个导数。 二阶导数在每个区间内线性变化。 所以,例如,如果你有

 xy y2 0 10 -30 2 5 -15 4 -5 -10 

那么x=1二阶导数是y2=-22.5 ,它在-30和-15之间。

你也可以创建一个新的splint()函数来返回二阶导数a*y2a[i]+b*y2a[i+1]