查看“數值微分”的源代码

{{NoteTA|G1=Math}}
'''數值微分'''是[[數值方法]]中的名詞，是用函數的值及其他已知資訊來估計一[[函數]][[導數]]的[[演算法]]。

[[Image:Derivative.svg|230px|right]]

==有限差分法==
最簡單的方式是使用[[有限差分]]近似。

簡單的二點估計法是計算經過(''x'',''f(x)'')及鄰近點(''x+h'',''f(x+h)'')二點形成[[割線]]的斜率<ref>Richard L. Burden, J. Douglas Faires (2000), ''Numerical Analysis'', (7th Ed),  Brooks/Cole. ISBN 0-534-38216-9</ref>選擇一個小的數值''h''，表示''x''的小變化，可以是正值或是負值。其斜率為
:<math>{f(x+h)-f(x)\over h}.</math>
此表示法是[[艾萨克·牛顿|牛頓]]的[[差商]]，也稱為一階[[均差]]。

割線斜率和切線斜率有些差異，差異大約和''h''成正比。若''h''近似於0，則割線斜率近似於切線斜率。因此，函數'''f'''真正在''x''處真正的斜率是割線趨近切線時的差商：
:<math>f'(x)=\lim_{h\to 0}{f(x+h)-f(x)\over h}.</math>

若直接將''h''用0取代會得到[[除以零]]的結果，因此計算導數需要一些較不直覺的的方式。

同様的，切線斜率也可以用(x - h)和x二點的割線斜率近似。

另外一種二點估計法是用經過(''x-h'',''f(x-h)'')和(''x+h'',''f(x+h)'')二點的割線，其斜率為
:<math>{f(x+h)-f(x-h)\over 2h}.</math>

上述公式稱為[[對稱差分]]，其一次項誤差相消，因此割線斜率和切線斜率的差和<math>h^2</math>成正比。對於很小的''h''而言這個值比單邊近似還要準確。特別的是公式雖計算x點的斜率，但不會用到函數在x點的數值。

估計誤差為：

:<math>R = {{-f^{(3)}(c)}\over {6}}h^2</math>,

其中<math>c</math>為在<math>x-h</math>和<math>x+h</math>之間的某一點。此誤差沒有包括因為有限準確度而產生的[[捨入誤差]]。

很多工程計算機都是用對稱差分來計算導數，像德州儀器（TI）的[[TI-82]]、[[TI-83]]、[[TI-84]]及[[TI-85]]，其''h''=0.001<ref name="Merseth2003">{{cite book|author=Katherine Klippert Merseth|title=Windows on Teaching Math: Cases of Middle and Secondary Classrooms|year=2003|publisher=Teachers College Press|isbn=978-0-8077-4279-2|page=34}}</ref><ref name="RubySellers2014">{{cite book|author1=Tamara Lefcourt Ruby|author2=James Sellers|author3=Lisa Korf |author4=Jeremy Van Horn |author5=Mike Munn|title=Kaplan AP Calculus AB & BC 2015|year=2014|publisher=Kaplan Publishing|isbn=978-1-61865-686-5|page=299}}</ref>。

雖然在實務十分常用，但上述二種方式的數值微分常被研究者批評，尤其是被一些鼓勵使用[[自動微分]]的研究者批評<ref name="GriewankWalther2008">{{cite book|author1=Andreas Griewank|author2=Andrea Walther|title=Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Second Edition|url=http://books.google.com/books?id=qMLUIsgCwvUC&pg=PA2|year=2008|publisher=SIAM|isbn=978-0-89871-659-7|pages=2–}}</ref>，因為上述的數值微分其精確度不高，若計算器精準度是六位數，用對稱差分計算導數只有三位數的精確度，而若是找到一計算斜率的函數，仍可以有幾乎六位數的精確度。例如假設f(x) = x<sup>2</sup>，用2x計算斜率有幾乎完整的準確度，而用差分近似就會有上述的問題。

===利用浮點數的實際考量===
[[Image:AbsoluteErrorNumericalDifferentiationExample.png|thumb|300px|在浮點數運算下，不同的<math>h</math>造成的捨入誤差及公式誤差，只有在特定值下誤差才是最小值]]

若計算時使用[[浮點數]]，就需要考慮''h''要取到多小。若選的太小，相減之後會有大的[[捨入誤差]]，事實上整個有限差分的公式都是[[条件数|病態]]的<ref name=Fornberg1>Numerical Differentiation of Analytic Functions, B Fornberg - ACM Transactions on Mathematical Software (TOMS), 1981</ref>，若''h''夠小，導數不為零的情形下，在相消後會得到數值微分為零的結果<ref name=SquireTrapp1>Using Complex Variables to Estimate Derivatives of Real Functions, W Squire, G Trapp - SIAM REVIEW, 1998</ref>，若''h''太大，計算割線斜率的結果就會更加準確，但用割線斜率估算切線斜率的誤差就更大了。

一種可以產生夠小的''h''，但又不會產生捨入誤差的方式是<math>\sqrt{\varepsilon}x</math>（不過''x''不能為0），其中{{le|最小浮点数|machine epsilon}}''&epsilon;''大約是2.2&times;10<sup>−16</sup>數量級。
<ref>Following ''Numerical Recipes in C'', [http://www.nrbook.com/a/bookcpdf/c5-7.pdf Chapter 5.7]</ref>。以下是一個一個可以平衡捨入誤差和公式誤差，有最佳精確度的''h''為

<math>h = 2\sqrt{\varepsilon|{f(x)\over f''(x)}|}</math>
<ref>p. 263 [http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h10/kompendiet/kap11.pdf]</ref>
（不過f"(x) = 0時不成立），而且需要有關函數的資訊。

上述的最小浮点数是針對雙精度（64-bit）變數，單精度變數在這類計算幾乎不太實用。其計算結果在二進制中不太可能是「整數」。雖然''x''是可以用浮點數表示的數字，但''x'' + ''h''幾乎不會也是可用浮點數表示（而且和''x''不同）的數字，因此''x'' + ''h''需調整為機器可讀的數字，因此會出現(''x'' + ''h'') - ''x'''''不等於'''''h''的情形，因此用二個函數計算值計算微分時，二個位置的差不會是''h''。幾乎所有的十進制分數在二進制下都會是循環小數（都像1/3在十進制中的情形一様），例如''h'' = 0.1在二進制下會是循環小數，是 0.000110011001100...。因此在浮點數下一個可能計算的方式是：
  h:=sqrt(eps)*x;
  xph:=x + h;
  dx:=xph - x;
  slope:=(F(xph) - F(x))/dx;

先計算(''x'' + ''h'') - ''x''的值，再用這個值作為微分算式的分母，不過若是用電腦計算，{{le|編譯器最佳化|Optimizing compiler}}的機能可能會認為''dx''和''h''相同，因此讓上述的方式失效。若是用C或其他類似的程式語言，可以讓''xph''宣告成[[Volatile变量]]，以避免此一問題。

===高階方法===
<!--{{further|Finite difference coefficients}}-->
也有用更高階估計導數的方法，或是估計高階導數的方法。

以下就是一階導數的五點法（一維下的{{le|五點模版|five-point stencil}}）<ref>Abramowitz & Stegun, Table 25.2</ref>
:<math>f'(x) = \frac{-f(x+2 h)+8 f(x+h)-8 f(x-h)+f(x-2h)}{12 h}+\frac{h^4}{30}f^{(5)}(c)</math>
其中<math>c\in[x-2h,x+2h]</math>.

==微分求积==
{{le|微分求积|Differential quadrature}}（Differential quadrature）是用函數在特定位置數值的加權和來近似導數<ref>Differential Quadrature and Its Application in Engineering: Engineering Applications, Chang Shu, Springer, 2000, ISBN 978-1-85233-209-9</ref><ref>Advanced Differential Quadrature Methods, Yingyan Zhang, CRC Press, 2009, ISBN 978-1-4200-8248-7</ref>，其名稱類似[[數值積分]]中用的求积（quadrature），也就是像[[梯形法]]或是[[辛普森法]]中用的加權和，有許多方式可找出加權的係數，在求解[[偏微分方程]]時會用到微分求积。

==複變的方法==
傳統用有限差分近似數值微分的方式是病態的，不過若<math>f</math>是[[全純函數]]，在實軸上的值都是實數，可以用複數平面中靠近<math>x</math>的位置來求值，此方式為[[數值穩定]]的方式，例如<ref name=SquireTrapp1/>一階導數可以用以下的複數導數公式計算<ref>{{cite journal | last1 = Martins | first1 = JRRA | first2 = P | last2 = Sturdza | first3 = JJ | last3 = Alonso | year = 2003 | id = {{citeseerx|10.1.1.141.8002}} | title = The Complex-Step Derivative Approximation | journal = ACM Transactions on Mathematical Software | volume = 29 | issue = 3 | pages = 245–262 | doi=10.1145/838250.838251}}</ref>：

:<math>f'(x)\approx \Im(f(x + ih))/h</math>.

上述公式只在計一階導數時有效，若要拓展到任意階導數，需要用到[[多重复数]]，結果也會是多重复数的導數。<ref>{{cite web |url=http://russell.ae.utexas.edu/FinalPublications/ConferencePapers/2010Feb_SanDiego_AAS-10-218_mulicomplex.pdf |title=存档副本 |accessdate=2012-11-24 |deadurl=yes |archiveurl=https://web.archive.org/web/20140109145840/http://russell.ae.utexas.edu/FinalPublications/ConferencePapers/2010Feb_SanDiego_AAS-10-218_mulicomplex.pdf |archivedate=2014-01-09 }}</ref>

而任意階的導數可以用[[柯西積分公式]]計算：
:<math>f^{(n)}(a) = {n! \over 2\pi i} \oint_\gamma {f(z) \over (z-a)^{n+1}}\, \mathrm{d}z</math>,
其中積分會用[[數值積分]]計算。

Lyness和Moler在1967年提出用複變數來計算數值微分<ref name=LynessMoler1>{{cite journal | first1 = J. N. | last1 = Lyness | first2 = C. B. | last2 = Moler | title = Numerical differentiation of analytic functions | journal = SIAM J.Numer. Anal. | volume = 4 | year = 1967 | pages = 202–210 | doi=10.1137/0704019}}</ref>。Abate和Dubner提出一種用複數[[拉普拉斯轉換]]的數值反演為基礎的算法<ref>{{cite journal | title = A New Method for Generating Power Series Expansions of Functions | first1 = J | last1 = Abate | first2 = H | last2 = Dubner | journal = SIAM J. Numer. Anal. | volume =5 | issue = 1 | pages = 102–112 |date=March 1968 | doi=10.1137/0705008}}</ref>。<!--An algorithm which can be used without requiring knowledge about the method or the character of the function was developed by Fornberg.<ref name=Fornberg1/>-->

==參考資料==
{{reflist}}

==相關條目==
*[[自動微分]]
*[[差分]]
*{{le|五点模板|Five-point stencil}}
*[[數值積分]]
*{{le|數值常微分方程|Numerical ordinary differential equations}}
*{{le|Savitzky–Golay濾波器|Savitzky–Golay filter}}
*{{le|數值分析軟體列表|List of numerical analysis software}}

== 外部連結 ==
* http://mathworld.wolfram.com/NumericalDifferentiation.html
* https://web.archive.org/web/20130820223117/http://math.fullerton.edu/mathews/n2003/NumericalDiffMod.html
*[http://numericalmethods.eng.usf.edu/topics/continuous_02dif.html Numerical Differentiation Resources: Textbook notes, PPT, Worksheets, Audiovisual YouTube Lectures] at [https://web.archive.org/web/20060906070428/http://numericalmethods.eng.usf.edu/ Numerical Methods for STEM Undergraduate]
*ftp://math.nist.gov/pub/repository/diff/src/DIFF Fortran code for the numerical differentiation of a function using Neville's process to extrapolate from a sequence of simple polynomial approximations.
* [http://www.nag.co.uk/numeric/fl/nagdoc_fl24/html/D04/d04conts.html NAG Library numerical differentiation routines]
* http://graphulator.com [http://graphulator.com Online numerical graphing calculator with calculus function.]

[[Category:数值分析]]
[[Category:微分学]]