复阶微分是一种利用复算法求实变量的实解析函数的一阶导数的数值的方法，避免了传统有限差分法固有的精度损失。

内容

刺激

这篇文章是在几个人的推动下写出来的。

几个月前，我在MathWorks餐厅遇到了Ricardo Pax Paxson，他是我们计算生物学小组的负责人。

然后，去年2月我在暹罗CS e会议期间，密歇根州大学航空航天工程教授Joaquim Martins。他给了我三篇论文，他写了关于他和他的同事使用的复杂步骤方法，以计算耦合计算流体动力学对飞机设计应用的结构有限元分析的模拟衍生物

Joaquim的论文引用了James Lyness和我发表在数值分析学报在1967年。那差不多是50年前的事了。

莱茵斯和硅藻土

詹姆斯·莱尼斯是我的朋友。我们是在苏黎世联邦理工学院认识的，当时我在那里做博士后。他是一个比我大几岁的英国人，大约在同一时间访问了ETH。那次访问之后，他在阿贡实验室度过了余生的职业生涯。他没有在LINPACK和EISPACK上工作过，但当我去Argonne工作时，我会在那里见到他。我经常展示的那张LINPACK的四位作者懒洋洋地躺在Jack Dongarra的车后座上的照片是在Lyness家的一次聚会上拍摄的。

那一年在苏黎世的研究成果非常丰富。詹姆斯有很多有趣的想法，我们发表了几篇论文。其中之一是《解析函数的数值微分》。我想这是第一次，介绍了实变量函数的复算方法。但实际的算法相当奇特。它们是由柯西积分定理推导出来的包括莫比乌斯数和欧拉麦克劳林求和公式。如果你今天要我详细描述它们，我是说不出来的。

我将要描述的复杂步骤算法比1967年那篇论文中描述的要简单得多。所以，尽管我们有一些优先权在实函数上使用复算的想法，我们当然不能说发明了今天的算法。

复杂的步骤算法首先是威廉郡和乔治特拉普的简介中描述的暹罗审查1998年。马丁与他的同事和学生一起提供了一个在下面的论文中的广泛指导。

该算法

我们关心的是分析函数。数学上，这意味着函数是无限可微的并且可以平滑地扩展到复平面上。在计算上，它可能意味着它是由一个“一行”公式定义的，而不是一个更广泛的代码段如果语句和为循环。

让$ f（z）$是这样的函数，让$ x_0 $是真实轴上的一个点，让$ h $是一个真实参数。将$ F（Z）$展开真实轴的泰勒序列。

$ $ (x_0 + ih) = F(从)+ ihF(从)- h ^ 2 F”(从)/ 2 ! ih ^ 3 F ^{(3)} / 3 ! +……＄$

采取双方的虚部，并将其分开以$ h $。

$ F'(x_0) = m(F(x_0+ih))/h + O(h^2) $

简单地计算函数$F$在虚参数$x_0+ih$处的值，然后除以$h$，就得到了导数$F'(x_0)$的近似值，精确到$O(h²)$。我们也可以选择$h=10^{-8}$。那么，近似中的误差大约与存储双精度浮点值$F'(x_0)$所涉及的舍入误差的大小相同。

这里是复杂的步长微分算法的整体:

$ h = 10^{-8} $

$ F'(x_0) \approx Im(F(x_0+ih))/h $$

一个例子

这是我最喜欢的例子之一。我在1967年纸上用过它，它已成为这项业务的一种标准。

$$ f（x）= \ frac {\ mbox {e} ^ {x}} {（\ cos {x}）^ 3 +（\ sin {x}）^ 3} $$

F = @(x) exp(x)./((cos(x))。^ 3 + (sin (x)) ^ 3)。

F = @ (x) exp (x) / ((cos (x)) ^ 3 + (sin (x)) ^ 3)。

这是$ x $ -axis的一部分的剧情，带有$ \ pi / 4 $的点。

ezplot (F[-π/ 4,π/ 2])轴([-π/ 4,π/ 2 0 6])集(gca,“xtick”[-π/ 4 0π/ 4π/ 2))((π/ 4π/ 4],[(π/ 4),F(π/ 4)),“标记”，“。”，'Markersize', 18)

求点处的斜率，也就是$F'(\ /4)$。

复杂的步长近似我还没有确定步长h是

Fpc = @(x,h) imag(F(x+i*h))/h;

步长中心有限差分近似h是

Fpd = @(x,h) (F(x+h) - F(x-h))/(2*h);

比较它们的递减序列h．

格式长disp ([h复阶有限差分]）为h = 10 ^ (- (1:16)) disp ([h Fpc(π/ 4,h) Fpd(π/ 4,h)))结束

H络合物步骤有限differerence 0.100000000000000 3.144276040634560 3.061511866568119 0.010000000000000 3.102180075411270 3.101352937655877 0.001000000000000 3.101770529535847 3.101762258158169 0.000100000000000 3.101766435192940 3.101766352480162 0.000010000000000 3.101766394249620 3.101766393398542 0.000001000000000 3.101766393840188 3.101766393509564 0.000000100000000 3.101766393836091 3.101766394841832 0.000000010000000 3.101766393836053 3.101766443691645 0.000000001000000 3.101766393836052 3.101766399282724 0.000000000100000 3.101766393836052 3.101763290658255 0.000000000010000 3.101766393836053 3.101785495118747 0.000000000001000 3.101766393836053 3.101963130802687 0.000000000000100 3.101766393836052 3.097522238704187 0.000000000000010 3.101766393836053 3.108624468950438 0.000000000000001 3.101766393836052 3.1086244689504380.00000000000000000 3.101766393836053 8.881784197001252

如预测所示，复杂的步骤已经通过序列获得了一半的全双精度结果H.．另一方面，有限差分法从来没有达到一半以上的精度加倍精度，完全崩溃的时间h达到最后一个值。

象征性的工具箱

由于我们有权访问符号工具箱，我们可以获得确切的答案。重新定义x和F成为象征性的。

信谊xF = exp(x)/(cos(x))^3 + (sin(x))^3

f = exp（x）/（cos（x）^ 3 + sin（x）^ 3）

求导并求在$\ /4$处的导数。

例如:Fps = simple(diff(F)) exact = subs(Fps,pi/4) flexact = double(exact)

Fps = (exp (x) *(因为罪(3 * x) + (3 * x) / 2 + (3 * sin (x)) / 2)) / (cos (x) ^ 3 + sin (x) ^ 3) ^ 2的= 2 ^ (1/2)* exp(π/ 4)flexact = 3.101766393836051

这验证了我们通过复步得到的结果。

错误的阴谋

现在我们得到了准确的答案，我们可以用复步和有限差分来计算一阶导数的误差。

h = 10 ^(-(1:16)。');犯错误= 0(16日2);为k = 1:16 errs(k,1) = abs(Fpc(pi/4,h(k))-exact);错(k, 2) = abs (Fpd(π/ 4 h (k))准确);结束logglog (h,errs)轴([eps 1 eps 1])集合(gca，“xdir”，“牧师”)传说(的复杂的步骤，'有限差异'，“位置”，“西南”)包含('梯级h') ylabel ('错误'）

这证实了复杂步长对于小步长是精确的，而有限差分方法永远无法达到完全的精度。

参考

《解析函数的数值微分》，4 .数值分析， 1967，页202-210。epubs.siam.org/doi/abs/10.1137/0704019

威廉·斯夸尔，乔治·特拉普，《用复变量估计实函数的导数》暹罗评论40.，1998年，第110-112页。epubs.siam.org/doi/abs/10.1137/S003614459631241X

《复阶导数近似法》，美国计算机学会数学软件汇刊， 2003, pp. 245-262，Martins2003CSD.pdf

Martins，Joaquim R. A.，复杂阶梯衍生物近似的指南，复杂的步骤

Kenway，Gaetan K.W.，Kennedy，Graeme J.和Martins，Joaquim R. A.，一种可扩展的高保真稳态空气弹性分析和伴随衍生计算，AIAA Journal（印刷机）

获取MATLAB代码

发布与MATLAB®R2013b

发布与MATLAB®R2013b