java 求导函数_在MATLAB中计算数值导数的最佳方法是什么?

这些只是一些快速而肮脏的建议 . 希望有人会发现它们有用!

1. Do you have a symbolic function or a set of points?

如果您有符号功能,您可以分析计算衍生物 . (很可能,如果那么容易的话,你会做到这一点,你不会在这里寻找替代品 . )

如果您具有符号函数且无法通过分析计算导数,则可以始终在一组点上评估函数,并使用此页面上列出的其他方法来评估导数 .

在大多数情况下,你有一组点(xi,fi),并且必须使用以下方法之一....

2. Is your grid evenly or unevenly spaced?

如果网格间隔均匀,您可能需要使用有限差分格式(请参阅维基百科文章here或here),除非您使用的是周期性边界条件(见下文) . Here是在解决网格上的常微分方程的背景下对有限差分方法的一个不错的介绍(特别参见幻灯片9-14) . 这些方法通常在计算上有效,易于实现,并且该方法的误差可以简单地估计为用于导出它的泰勒展开的截断误差 .

如果您的网格间距不均匀,您仍然可以使用有限差分格式,但表达式更加困难,并且精度随着网格的统一程度而变化很大 . 如果您的网格非常不均匀,您可能需要使用较大的模板尺寸(更多相邻点)来计算给定的导数点 . 人们经常构造插值多项式(通常是Lagrange polynomial)并区分该多项式来计算导数 . 例如,参见this StackExchange问题 . 使用这些方法通常很难估计错误(尽管有些人试图这样做:here和here) . Fornberg's method在这些情况下通常非常有用....

必须注意您域名的边界,因为模板通常涉及域外的点 . 有些人引入"ghost points"或将边界条件与不同阶数的导数组合以消除这些并简化模板 . 另一种方法是使用右侧或左侧有限差分方法 .

Here's一种优秀的有限差分方法,包括低阶的中心,右侧和左侧方案 . 我在工作站附近打印出这个,因为我发现它非常有用 .

3. Is your domain periodic? Can you assume periodic boundary conditions?

如果您的域是周期性的,您可以使用傅里叶谱方法计算导数到非常高的阶数 . 这种技术在某种程度上牺牲了性能以获得高精度 . 实际上,如果您使用的是N点,那么您对导数的估计大约是第N阶的准确度 . 有关更多信息,请参阅(例如)this WikiBook .

傅立叶方法通常使用快速傅立叶变换(FFT)算法来实现大致O(N log(N))性能,而不是天然实现的离散傅立叶变换(DFT)可能采用的O(N ^ 2)算法 .

如果您的函数和域不是周期性的,则不应使用傅里叶谱方法 . 如果您尝试将其与非周期性函数一起使用,则会出现大的错误和不良的"ringing"现象 .

计算任何阶数的导数需要1)从网格空间到光谱空间的变换(O(N log(N))),2)傅里叶系数乘以其光谱波数(O(N)),和2)从光谱空间到网格空间的逆变换(再次为O(N log(N))) .

在将傅里叶系数乘以其光谱波数时必须小心 . FFT算法的每个实现似乎都有自己的频谱模式和归一化参数的排序 . 例如,请参阅Math StackExchange上this question的答案,以获取有关在MATLAB中执行此操作的说明 .

4. What level of accuracy are you looking for? Do you need to compute the derivatives within a given tolerance?

出于许多目的,一阶或二阶有限差分方案可能就足够了 . 为了获得更高的精度,您可以使用更高阶的泰勒展开式,从而降低高阶项 .

如果需要计算给定容差范围内的导数,可能需要查看具有所需错误的高阶方案 .

通常,减少误差的最佳方法是在有限差分方案中减少网格间距,但这并非总是可行的 .

请注意,高阶有限差分格式几乎总是需要更大的模板尺寸(更多的相邻点) . 这可能会导致边界问题 . (参见上面关于鬼点的讨论 . )

5. Does it matter to you that your derivative is evaluated on the same points as your function is defined?

MATLAB提供 diff 函数来计算相邻数组元素之间的差异 . 这可以用于通过一阶前向差分(或前向有限差分)方案来计算近似导数,但估计是低阶估计 . 如MATLAB的 diff (link)文档中所述,如果输入长度为N的数组,它将返回一个长度为N-1的数组 . 当您在N点上使用此方法估计导数时,您只能估算出N-1点的导数 . (请注意,如果它们按升序排序,则可以在不均匀网格上使用 . )

在大多数情况下,我们希望在所有点评估衍生物,这意味着我们想要使用除 diff 方法之外的其他东西 .

6. Do you need to calculate multiple orders of derivatives?

可以 Build 一个方程组,其中网格点函数值和这些点上的一阶和二阶导数都相互依赖 . 这可以通过像往常一样在相邻点处组合泰勒展开来找到,但是保持导数项而不是将它们抵消,并将它们与相邻点的那些连接在一起 . 这些方程式可以通过线性代数求解,不仅可以得到一阶导数,也可以得到第二个导数(或者更高阶,如果设置得当) . 我相信这些被称为组合有限差分方案,它们通常与紧致有限差分方案结合使用,这将在下面讨论 .

紧凑有限差分格式(link) . 在这些方案中, Build 一个设计矩阵并通过矩阵求解同时计算所有点的导数 . 它们被称为"compact"因为它们通常被设计为需要比具有可比精度的普通有限差分方案更少的模板点 . 因为它们涉及将所有点连接在一起的矩阵方程,确定据说紧凑有限差分方案具有"spectral-like resolution"(例如Lele's 1992 paper - 优秀!),这意味着它们通过依赖于所有节点值来模拟频谱方案,并且因此,它们在所有长度尺度上保持准确性 . 相反,典型的有限差分方法仅局部精确(例如,点#13处的导数通常不依赖于点#200处的函数值) .

目前的研究领域是如何在紧凑的模板中最好地解决多种衍生物 . 这种研究,组合,紧凑的有限差分方法的结果是强大的和广泛适用的,尽管许多研究人员倾向于针对特定需求(性能,准确性,稳定性或诸如流体动力学的特定研究领域)调整它们 .

Ready-to-Go Routines

如上所述,可以使用 diff 函数(文档link)来计算相邻数组元素之间的粗略导数 .

MATLAB的 gradient 例程(link到文档)是一个很好的选择,用于许多目的 . 它实现了二阶中心差分方案 . 它具有计算多维导数并支持任意网格间距的优点 . (感谢@thewaywewalk指出这个明显的遗漏!)

我用Fornberg的方法(见上文)开发了一个小程序( nderiv_fornberg )来计算任意网格间距的一维的有限差分 . 我觉得它很容易使用 . 它在边界处使用6个点的侧面模板,在内部使用中心的5点模板 . 它可以在MATLAB File Exchange here中找到 .

Conclusion

数值微分领域非常多样化 . 对于上面列出的每种方法,有许多变体具有各自的优点和缺点 . 这篇文章几乎不是数值微分的完整处理 .

每个应用程序都不同 . 希望这篇文章能够为感兴趣的读者提供有组织的考虑因素和资源清单,以便选择适合自己需要的方法 .

可以使用特定于MATLAB的代码片段和示例来改进此社区wiki .

MATLAB处理不等距数值点进行导数计算,可以通过多种方式实现。一种常见的方法是使用数值微分公式,比如有限差分法。有限差分法可以是前向差分、后向差分或心差分。对于不等距数据点,通常需要采用加权差分法或者使用插值后再计算导数方法。 一种方法是使用插值来估计连续函数,然后再在这些插值点上计算导数。例如,可以使用MATLAB的`interp1`函数进行插值,并结合差分公式来计算导数。`interp1`函数可以根据不等距的x值,返回等距的y值,然后可以在这些等距的y值上使用`diff`函数来近似导数。 另外,MATLAB的符号计算工具箱(Symbolic Math Toolbox)提供了一些函数,可以直接在不等距数据上计算导数。例如,`gradient`函数可以接受一个函数表达式和一组点作为输入,然后计算出在这些点上的导数值。对于不等距的数据点,可以将x值和y值作为向量传给`gradient`函数,从而得到导数的近似值。 下面是一个简单的例子: ```matlab % 假设x和y是已知的不等距数值点向量 x = [x1, x2, x3, ..., xn]; y = [y1, y2, y3, ..., yn]; % 对y进行插值,这里使用线性插值作为例子 yi = interp1(x, y, linspace(min(x), max(x), 1000), 'linear'); % 计算插值点上的导数 % 注意:这里假设x和y是等距的,因为我们使用了linspace生成等距的插值点 dyi = gradient(yi, linspace(min(x), max(x), 1000)); % 如果需要对特定的不等距点计算导数,可以考虑使用数值微分公式 % 如心差分公式: dx = diff(x); dy = diff(y); % 心差分近似导数 dydx = (dy(1:end-1) + dy(2:end)) ./ (dx(1:end-1) + dx(2:end)); ``` 这种方法结合了插值和差分来处理不等距数据点的导数计算。需要注意的是,导数的准确性受到插值方法和数据点分布的影响。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值