如何计算在无穷远处涉及指数和奇点的函数的数值二阶导数。 不幸的是,Ridder的方法在“C中的数值配方”中提供的数值导数只能计算一阶导数(它需要预先对函数进行解析表达式。)此外,我尝试了Chebyshev近似并在之后区分函数,但给出的值是方式关闭实际值。 我也尝试过在数学论文中提供的一些有限差分算法,但它们也容易出错。 函数是e ^(x / 2)/ x ^ 2。 对此事我感激不尽。
提前致谢
最新编辑:问题解决了C ++中提供的FADBAD库做得非常好。 它们可以通过https://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, "nn"); 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, "nn"); 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)。 它会在任何位置为您提供第二个衍生物。
使用x
和y
值调用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]
以上就是c/c++开发分享数值微分相关内容,想了解更多C/C++开发(异常处理)及C/C++游戏开发关注计算机技术网(www.ctvol.com)!)。
本文来自网络收集,不代表计算机技术网立场,如涉及侵权请联系管理员删除。
ctvol管理联系方式QQ:251552304
本文章地址:https://www.ctvol.com/c-cdevelopment/559843.html