如何在matlab中评估函数的导数？

Question

这应该很简单.我有一个函数f(x),我想在MATLAB中f'(x)给出一个给定x的函数.

我的所有搜索都提出了符号数学,这不是我需要的,我需要数字微分.

例如,如果我定义: fx = inline('x.^2')

我想找到说f'(3),这是6,我不想找2x

Answer 1

如果已知您的函数是两次可微分,请使用

f'(x) = (f(x + h) - f(x - h)) / 2h

这是h的二阶精度.如果它只是一次可区分,请使用

f'(x) = (f(x + h) - f(x)) / h     (*)

这是h的第一顺序.

这是理论.在实践中,事情非常棘手.我将采用第二个公式(第一个顺序),因为分析更简单.将第二个订单作为练习.

第一个观察是你必须确保(x + h) - x = h,否则你会得到巨大的错误.实际上,f(x + h)和f(x)彼此接近(比如2.0456和2.0467),当你减去它们时,你会失去很多重要的数字(这里是0.0011,有3位有效数字少)比x).因此,h上的任何错误都可能对结果产生巨大影响.

所以,第一步,修复一个候选人h(我会在一分钟内告诉你如何选择它),并将你的计算数量h'=(x + h) - x.如果您使用的是C语言,则必须注意将h或x定义为易失性,以便不进行优化.

接下来,选择h.(*)中的错误有两部分:截断错误和舍入错误.截断错误是因为公式不精确:

(f(x + h) - f(x)) / h = f'(x) + e1(h)

哪里e1(h) = h / 2 * sup_{x in [0,h]} |f''(x)|.

舍入误差来自f(x + h)和f(x)彼此接近的事实.它可以大致估计为

e2(h) ~ epsilon_f |f(x) / h|

其中epsilon_ff(x)(或f(x + h)的计算中的相对精度是近的).这必须从您的问题进行评估.对于简单的功能,epsilon_f可以作为机器epsilon.对于更复杂的,它可能比数量级更糟.

所以你想要h最小化e1(h) + e2(h).将所有内容整合在一起并优化h产量

h ~ sqrt(2 * epsilon_f * f / f'')

这必须从你的功能估计.你可以粗略估计一下.如有疑问,请使用h~sqrt(epsilon),其中epsilon =机器精度.对于h的最佳选择,导数已知的相对精度是sqrt(epsilon_f),即.有效数字的一半是正确的.

简而言之:太小啊=>舍入错误,太大啊=>截断错误.

对于二阶公式,相同的计算产生

h ~ (6 * epsilon_f / f''')^(1/3)

并且导数的(epsilon_f)^(2/3)的分数精度(假设双精度,它通常比一阶公式好一个或两个有效数字).

如果这太不精确,可以随意提出更多方法,有很多技巧可以提高准确性.Richardson推断是平滑函数的良好开端.但是这些方法通常会计算很多次,如果你的函数很复杂,这可能或不是你想要的.

如果你要在不同的点上多次使用数值导数,那么构造切比雪夫近似就变得很有趣了.

Answer 2

要获得数值差异(对称差异),您需要计算 (f(x+dx)-f(x-dx))/(2*dx)

fx = @(x)x.^2;
fPrimeAt3 = (fx(3.1)-fx(2.9))/0.2;

或者,您可以创建函数值向量并应用DIFF,即

xValues = 2:0.1:4;
fValues = fx(xValues);
df = diff(fValues)./0.1; 

请注意,它diff采用前向差异,并假设dx等于1.

但是,在您的情况下,最好定义fx为多项式,并评估函数的导数,而不是函数值.

Answer 3

缺乏符号工具箱,没有什么能阻止你使用Derivest,这是一种自动自适应数值微分的工具.

derivest(@sin,pi)
ans =
           -1

对于你的例子,它非常好.实际上,它甚至可以估算得到的近似误差.

fx = inline('x.^2');
[fp,errest] = derivest(fx,3)

fp =
            6
errest =
   3.6308e-14