在MATLAB中計算數值導數的最佳方法是什麼？

Question

（注：這是爲了成爲一個社區維基）

假設我有一組點喜 = {X0，X1，X2 ... XN }和相應的函數值音響 = F（XI）= {F0，F1，F2，...，FN }，其中f（x）通常是未知函數。（在某些情況下，我們可能知道˚F（X）時間提前，但我們要普遍做到這一點，因爲我們經常不知道˚F（提前X）。） 什麼是近似的衍生物的好方法？（x）在每個點xi？即，如何可以估算值DFI == d/d X音響 == d ˚F（XI）/ d X在每個點XI的？

不幸的是，MATLAB沒有很好的通用數值微分例程。部分原因可能是因爲選擇一個好的例程可能很困難！

那麼有什麼樣的方法？什麼樣的程序存在？我們如何爲特定問題選擇一個好的例程？

選擇如何在MATLAB區分時，有幾個方面的考慮：

1. 你有一個象徵性的功能或一組點？
2. 您的網格是否均勻或不均勻間隔？
3. 您的域名是否有周期性？你能假設週期性的邊界條件嗎？
4. 你在尋找什麼樣的準確性？你是否需要計算給定容差內的衍生物？
5. 對你而言，你的衍生物是否在你的函數被定義的相同點上進行評估對你有影響嗎？
6. 您是否需要計算衍生產品的多個訂單？

什麼是繼續進行的最佳方式？

Answer 1

這些只是一些快速和骯髒的建議。希望有人會發現他們有幫助！

1。你有一個象徵性的功能或一組點嗎？

• 如果你有一個象徵性的功能，你可以能夠分析計算出衍生物。 （很可能，如果這麼簡單，你就可以做到這一點，並且你不會在這裏尋找替代品。）
• 如果你有一個符號功能並且不能分析計算導數，你總是可以在點的集合，並使用此頁面上列出的其他方法來評估導數。
• 在大多數情況下，你有一組點（XI，FI），且必須使用下列方法之一....

2.您的網格均勻或不均勻間隔多長時間？

• 如果您的網格是均勻間隔，你可能會想用有限差分法（看到的要麼是維基百科的文章here或here的），除非你使用週期邊界條件（見下文）。 Here是在網格上求解常微分方程的上下文中的有限差分方法的體面介紹（特別參見幻燈片9-14）。這些方法通常計算效率高，實現簡單，並且該方法的誤差可以簡單地估計爲用於推導它的泰勒展開的截斷誤差。
• 如果你的網格是不均勻間隔，你仍然可以使用有限差分方案，但表達式更困難，精度因網格的均勻性而異。如果網格非常不均勻，則可能需要使用大模板大小（更多相鄰點）來計算給定點處的導數。人們經常構造一個內插多項式（通常是Lagrange polynomial）並區分該多項式來計算導數。例如，請參閱this StackExchange問​​題。使用這些方法估計錯誤通常很困難（儘管有些人試圖這樣做：here和here）。 Fornberg's method在這些情況下通常非常有用......
• 必須小心在您的域的邊界，因爲模板通常涉及域外的點。有些人會引入「鬼點」或將邊界條件與不同順序的衍生物相結合，以消除這些「鬼點」並簡化模板。另一種方法是使用右側或左側有限差分方法。
• Here's an 優秀有限差分法的「備忘單」，包括低階中心，右側和左側方案。我在我的工作站附近打印輸出，因爲我覺得它非常有用。

3.您的域名是否有周期性？你能假設週期性的邊界條件嗎？

• 如果您的域是週期性的，您可以使用傅里葉譜方法計算極高精度的導數。這種技術在某種程度上犧牲了性能以獲得高精度。事實上，如果你使用N個點，那麼你對導數的估計大約是N階精度。有關更多信息，請參閱（例如）this WikiBook。
• 傅里葉方法通常使用快速傅立葉變換（FFT）算法來實現大致O（N log（N））性能，而不是使用天真實現的離散傅立葉變換（DFT）可能會出現的O（N^2）算法採用。
• 如果你的函數和域是而不是週期性的，你不應該使用傅里葉譜方法。如果嘗試使用週期性的而不是的函數，將會出現大的錯誤和不希望出現的「振鈴」現象。計算任何次序的導數需要1）從網格空間到頻譜空間的變換（O（N log（N））），2）傅立葉係數乘以它們的頻譜波數（O（N）），以及2）從頻譜空間到網格空間的逆變換（再次O（N log（N）））。
• 將傅立葉係數乘以它們的光譜波數時必須小心。 FFT算法的每個實現似乎都有自己的頻譜模式和歸一化參數排序。例如，請參閱Math StackExchange上關於this question的回答，以獲取有關在MATLAB中執行此操作的說明。

4.你在尋找什麼準確度？你是否需要計算給定容差內的衍生物？

• 對於很多目的，第一或第二階有限差分方案可能就足夠了。爲了獲得更高的精度，可以使用更高階泰勒展開式，從而降低高階項。
• 如果您需要計算給定容差內的導數，則可能需要查找具有所需誤差的高階方案。
• 通常，減少錯誤的最佳方法是減少有限差分方案中的網格間距，但這並非總是可行。
• 請注意，高階有限差分方案几乎總需要更大的模板尺寸（更多相鄰點）。這可能會導致邊界問題。 （參見上述關於鬼點的討論。）

5.不要緊，你爲你的函數定義你的導數在同一個點評估？

• MATLAB提供了diff函數來計算相鄰數組元素之間的差異。這可用於通過一階正向差分（或正向有限差分）方案計算近似導數，但估計值是低階估計值。如MATLAB的文檔diff（link）所述，如果輸入長度爲N的數組，它將返回一個長度爲N-1的數組。當您在N個點上使用這種方法估計導數時，您只能估計N-1個點處的導數。 （請注意，這可以在不平坦的網格被使用，如果它們被按升序排序。）
• 在大多數情況下，我們希望在所有點評估的導數，這意味着我們要使用的東西除了diff方法。

6.您是否需要計算衍生產品的多個訂單？

• 可以建立一個方程組，其中網格點函數值和這些點上的一階和二階導數都相互依賴。這可以通過在鄰近點照常結合泰勒展開，但保持導數項，而不是取消他們，並與相鄰點的連接在一起，它們被發現。這些方程可以通過線性代數求解，不僅可以給出一階導數，還可以給出第二個（或者更高的階數，如果設置正確的話）。相信這些被稱爲組合有限差分格式，並且它們常常結合緊湊有限差分格式，這將在下面討論使用。
• 緊緻有限差分方案（link）。在這些方案中，建立一個設計矩陣並通過矩陣求解同時計算所有點的導數。它們被稱爲「緊湊型」，因爲它們通常被設計爲比相當準確度的普通有限差分方案需要更少的模板點。因爲它們涉及矩陣方程所有的點連接起來，某些緊湊有限差分法，據說有「譜像分辨率」（例如Lele's 1992 paper - ！優秀），這意味着它們將根據所有節點的值模擬頻譜方案並因此，它們在所有長度範圍內都保持準確。相反，典型的有限差分方法僅在局部精確（例如，點＃13的導數通常不取決於點＃200的函數值）。
• 研究當前地區是如何最好地解決在緊湊蠟紙多衍生物。儘管許多研究人員傾向於針對特定需求（性能，準確性，穩定性或諸如流體動力學的特定研究領域）調整它們，但這種研究的結果，組合，緊湊的有限差分方法是強大且廣泛適用的。

準備對去例程

• 如上所述，一個可以使用diff函數（link到文檔）來計算相鄰的陣列元件之間的粗衍生物。

• MATLAB的gradient程序（link到文檔）是用於多種用途的理想選擇。它實現了二階中心差分格式。它具有計算多維度導數和支持任意網格間隔的優點。 （感謝@thewaywewalk指出這個明顯的遺漏！）

• 我用Fornberg的方法（見上）開發了一個小例程（nderiv_fornberg）來計算任意網格間距的一維有限差分。我覺得它很容易使用。它在邊界處使用了6點的雙面模板，並在內部使用了一個居中的5點模板。它可以在MATLAB文件交換here。

結論

數值微分的領域是非常多樣的。對於上面列出的每種方法，有許多變體都有自己的一套優點和缺點。這篇文章幾乎不是數字分化的完整處理。

每個應用程序都是不同的。希望這篇文章能夠爲有興趣的讀者提供一個有組織的考慮因素和資源列表，供您選擇適合自己需求的方法。

此社區wiki可以通過代碼片段和MATLAB特有的示例進行改進。

Answer 2

我相信這些特定問題還有更多。所以我詳細闡述瞭如下主題：

（4）問：你在尋找什麼級別的準確度？你是否需要計算給定容差內的衍生物？

答：數值微分的準確性對於感興趣的應用是主觀的。通常，它的工作方式是，如果您在前向問題中使用ND來近似導數以根據感興趣的信號估計特徵，那麼您應該瞭解噪聲擾動。通常，這些僞影包含高頻分量，並且通過定義差分器，噪聲效應將以$ i \ omega^n $的幅度級被放大。所以，增加微分器的精度（增加多項式精度）根本無濟於事。在這種情況下，您應該能夠抵消噪音對差異化的影響。這可以按照匹配順序完成：首先使信號平滑，然後區分。但更好的方法是使用「低通微分器」。 MATLAB庫的一個很好的例子可以在here找到。但是，如果情況並非如此，並且您在反演問題（如solvign偏微分方程）中使用ND，則微分器的全局精度非常重要。根據什麼樣的彈性條件（BC）適合您的問題，設計將會相應地進行調整。重擊的規則是增加已知的全頻帶微分器的數值精度。您需要設計一個處理合適BC的衍生矩陣。您可以使用上述鏈接找到針對此類設計的全面解決方案。 （5）對你而言，你的衍生物是否在你的函數被定義的相同點上進行評估對你有影響嗎？答：絕對是。對相同網格點的ND進行評估稱爲「集中式」，而不是「交錯」方案。請注意，使用奇數階導數，集中ND會偏離差分器的頻率響應精度。因此，如果您在逆向問題中使用這種設計，這會干擾您的近似。而且，相反的情況適用於由交錯方案使用的偶數階微分的情況。您可以使用上面的鏈接找到關於此主題的全面解釋。

（6）您是否需要計算衍生產品的多個訂單？ 這完全取決於您的應用程序。您可以參考我提供的相同鏈接，並關注多個派生設計。