数值微分 [英] Numerical Differentiation
问题描述
如何计算涉及无穷远指数和奇点的函数的数值二阶导数.不幸的是，C 中的数值食谱"中提供的 Ridder 方法的数值导数只能计算一阶导数(它需要事先对函数进行解析表达式.)此外，我尝试过切比雪夫近似并在之后对函数进行微分，但给出的值偏离实际值.我还尝试了数学论文中提供的一些有限差分算法，但它们也很容易出错.函数是 e^(x/2)/x^2.我将不胜感激.
How can I calculate the numerical second derivative of a function involving an exponential and a singularity at infinity. Unfortunately, the numerical derivative by Ridder's methods provided in "Numerical Recipes in C" can only calculate the first derivative (It requires analytical expression of the function beforehand.) Furthermore I have tried Chebyshev approximation and differentiating the function afterwards but the values given were way off the actual values. I have also tried some finite difference algorithms provided in a mathematical paper yet they were error prone too. The function is e^(x/2) / x^2. I would appreciate any help on the matter.
提前致谢
最新问题已解决，C++ 中可用的 FADBAD 库做得非常好.它们可通过 http://www.fadbad.com/fadbad.html
Latest The issue was solved the FADBAD libraries available in C++ did an extremely good job. They are available via http://www.fadbad.com/fadbad.html
// The compilation command used is given below
// gcc Q3.c nrutil.c DFRIDR.c -lm -o Q3
#include <stdio.h>
#include <math.h>
#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 < b; i++)
{
x = a + i / pre;
y = func(x);
nd1 = dfridr(func, x, 0.1, &err);
//nd2 = dfridr((*func), x, 0.1, &err);
fprintf(fp, "%f \t %f \t %f \t %f \n", x, y, nd1);
if(y > 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 <stdio.h>
#include <math.h>
#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 < b; i++)
{
x = a + i / pre;
y = func(x);
nd1 = chebev(a,b,cder,NVAL,x);
nd2 = chebev(a,b,cder2,NVAL,x);
fprintf(fp, "%f \t %f \t %f \t %f \n", x, y, nd1, nd2);
if(y > 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;
}
推荐答案
该函数是可微的，因此使用数字方法可能不是最好的.二阶导数是:
That function is differentiable so using a numeric method is likely not the best. The second derivative is:
6*exp(x/2)/(x^4)-2*exp(x/2)/x^3 + exp(x/2)/(4*x^2)
6*exp(x/2)/(x^4)-2*exp(x/2)/x^3 + exp(x/2)/(4*x^2)
当然可以简化上面的内容以加快计算速度.第一次把原来的公式弄错了.
The above can be simplified of course to speed up computation. had original formula wrong the first time.
这篇关于数值微分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
