Mathematica生成高阶微分公式

本文包括常见的高阶的有限差分公式，以及使用 Mathematica 自动计算数值微分公式的代码。

原理就是基于待定系数法求解得到如下形式的有限差分公式

\[ f^{(\alpha)}(x) = \frac{1}{h^\alpha}\left( c_{-m} f(x-m h) + c_{-m+1} f(x-(m-1)h) +\cdots + c_n f(x+n h) \right) + O(h^\beta) \]

其中的\(m,n \ge 0\) 代表模板宽度，要求\(m+n \ge \alpha\)。

Mathematica 源代码

FDFormula[m_, n_, a_, {f_, x_, h_}] :=
 Module[{mat, b, coeff, expr, i, j},
  mat = Table[If[j == 0, 1, i^j], {j, 0, m + n}, {i, -m, n}];
  b = Table[If[i == a, a!, 0], {i, 0, n + m}];
  coeff = LinearSolve[mat, b];
  expr = Sum[coeff[[m + 1 + i]] f[x + i h ]/h^a, {i, -m, n}];
  For[j = a + 1, j <= m + n + 2, j++,
   If[Simplify@Sum[coeff[[m + 1 + i]]*i^j, {i, -m, n}] != 0, Break[]]];
  Print["coeff: ", coeff, "\nexpression: ",
   Format[expr, TraditionalForm], "\nseries: ",
   Series[expr, {h, 0, j - a}], "\norder: ", j - a];
  Print[Format[expr, TeXForm], "=",
   Format[Series[expr, {h, 0, j - a}], TeXForm]];
  Return[coeff]]

使用示例

例一，求解\(f'\)的三个点的中心差分公式

FDFormula[1, 1, 1, {f, x, h}]

这里的{f,x,h}要保持未定义的状态，返回结果为

coeff: {-(1/2),0,1/2}
expression: f(h+x)/(2 h)-f(x-h)/(2 h)
series: (f^\[Prime])[x]+1/6 (f^(3))[x] h^2+O[h]^3
order: 2

对应的是经典的二阶中心差分公式

\[ \frac{f(x+h)}{2 h}-\frac{f(x-h)}{2 h}=f'(x)+\frac{1}{6} h^2 f^{(3)}(x)+O\left(h^3\right) \]

例二，求解\(f'\)的五个点的中心差分公式

FDFormula[2, 2, 1, {f, x, h}]

这里的{f,x,h}要保持未定义的状态，返回结果为

coeff: {1/12,-(2/3),0,2/3,-(1/12)}
expression: f(x-2 h)/(12 h)-(2 f(x-h))/(3 h)+(2 f(h+x))/(3 h)-f(2 h+x)/(12 h)
series: (f^\[Prime])[x]-1/30 (f^(5))[x] h^4+O[h]^5
order: 4

对应的是四阶五点的中心差分公式

\[ \frac{f(x-2 h)}{12 h}-\frac{2 f(x-h)}{3 h}+\frac{2 f(x+h)}{3 h}-\frac{f(x+2 h)}{12 h}=f'(x)-\frac{1}{30} h^4 f^{(5)}(x)+O\left(h^5\right) \]

例三，求解\(f''\)的迎风差分公式

FDFormula[0, 4, 2, {f, x, h}]

这里的{f,x,h}要保持未定义的状态，返回结果为

coeff: {35/12,-(26/3),19/2,-(14/3),11/12}
expression: (35 f(x))/(12 h^2)-(26 f(h+x))/(3 h^2)+(19 f(2 h+x))/(2 h^2)-(14 f(3 h+x))/(3 h^2)+(11 f(4 h+x))/(12 h^2)
series: (f^\[Prime]\[Prime])[x]+5/6 (f^(5))[x] h^3+O[h]^4
order: 3

对应的迎风公式如下

\[ \frac{35 f(x)}{12 h^2}-\frac{26 f(x+h)}{3 h^2}+\frac{19 f(x+2 h)}{2 h^2}-\frac{14 f(x+3 h)}{3 h^2}+\frac{11 f(x+4 h)}{12 h^2}=f''(x)+\frac{5}{6} h^3 f^{(5)}(x)+O\left(h^4\right) \]

中心差分

下面提供可能会用到的若干个中心差分格式。

一阶导

三点（二阶）

\[ \frac{f(x+h)}{2 h}-\frac{f(x-h)}{2 h}=f'(x)+\frac{1}{6} h^2 f^{(3)}(x)+O\left(h^3\right) \]

五点（四阶）

\[ \frac{f(x-2 h)}{12 h}-\frac{2 f(x-h)}{3 h}+\frac{2 f(x+h)}{3 h}-\frac{f(x+2 h)}{12 h}=f'(x)-\frac{1}{30} h^4 f^{(5)}(x)+O\left(h^5\right) \]

七点（六阶）

\[ -\frac{f(x-3 h)}{60 h}+\frac{3 f(x-2 h)}{20 h}-\frac{3 f(x-h)}{4 h}+\frac{3 f(x+h)}{4 h}-\frac{3 f(x+2 h)}{20 h}+\frac{f(x+3 h)}{60 h}=f'(x)+\frac{1}{140} h^6 f^{(7)}(x)+O\left(h^7\right) \]

二阶导

三点（二阶）

\[ -\frac{2 f(x)}{h^2}+\frac{f(x-h)}{h^2}+\frac{f(x+h)}{h^2}=f''(x)+\frac{1}{12} h^2 f^{(4)}(x)+O\left(h^3\right) \]

五点（四阶）

\[ -\frac{5 f(x)}{2 h^2}-\frac{f(x-2 h)}{12 h^2}+\frac{4 f(x-h)}{3 h^2}+\frac{4 f(x+h)}{3 h^2}-\frac{f(x+2 h)}{12 h^2}=f''(x)-\frac{1}{90} h^4 f^{(6)}(x)+O\left(h^5\right) \]

七点（六阶）

\[ -\frac{49 f(x)}{18 h^2}+\frac{f(x-3 h)}{90 h^2}-\frac{3 f(x-2 h)}{20 h^2}+\frac{3 f(x-h)}{2 h^2}+\frac{3 f(x+h)}{2 h^2}-\frac{3 f(x+2 h)}{20 h^2}+\frac{f(x+3 h)}{90 h^2}=f''(x)+\frac{1}{560} h^6 f^{(8)}(x)+O\left(h^7\right) \]

三阶导

五点（二阶）

\[ -\frac{f(x-2 h)}{2 h^3}+\frac{f(x-h)}{h^3}-\frac{f(x+h)}{h^3}+\frac{f(x+2 h)}{2 h^3}=f^{(3)}(x)+\frac{1}{4} h^2 f^{(5)}(x)+O\left(h^3\right) \]

七点（四阶）

\[ \frac{f(x-3 h)}{8 h^3}-\frac{f(x-2 h)}{h^3}+\frac{13 f(x-h)}{8 h^3}-\frac{13 f(x+h)}{8 h^3}+\frac{f(x+2 h)}{h^3}-\frac{f(x+3 h)}{8 h^3}=f^{(3)}(x)-\frac{7}{120} h^4 f^{(7)}(x)+O\left(h^5\right) \]

九点（六阶）

\[ -\frac{7 f(x-4 h)}{240 h^3}+\frac{3 f(x-3 h)}{10 h^3}-\frac{169 f(x-2 h)}{120 h^3}+\frac{61 f(x-h)}{30 h^3}-\frac{61 f(x+h)}{30 h^3}+\frac{169 f(x+2 h)}{120 h^3}-\frac{3 f(x+3 h)}{10 h^3}+\frac{7 f(x+4 h)}{240 h^3}=f^{(3)}(x)+\frac{41 h^6 f^{(9)}(x)}{3024}+O\left(h^7\right) \]

四阶导

五点（二阶）

\[ \frac{6 f(x)}{h^4}+\frac{f(x-2 h)}{h^4}-\frac{4 f(x-h)}{h^4}-\frac{4 f(x+h)}{h^4}+\frac{f(x+2 h)}{h^4}=f^{(4)}(x)+\frac{1}{6} h^2 f^{(6)}(x)+O\left(h^3\right) \]

七点（四阶）

\[ \frac{28 f(x)}{3 h^4}-\frac{f(x-3 h)}{6 h^4}+\frac{2 f(x-2 h)}{h^4}-\frac{13 f(x-h)}{2 h^4}-\frac{13 f(x+h)}{2 h^4}+\frac{2 f(x+2 h)}{h^4}-\frac{f(x+3 h)}{6 h^4}=f^{(4)}(x)-\frac{7}{240} h^4 f^{(8)}(x)+O\left(h^5\right) \]

九点（六阶）

\[ \frac{91 f(x)}{8 h^4}+\frac{7 f(x-4 h)}{240 h^4}-\frac{2 f(x-3 h)}{5 h^4}+\frac{169 f(x-2 h)}{60 h^4}-\frac{122 f(x-h)}{15 h^4}-\frac{122 f(x+h)}{15 h^4}+\frac{169 f(x+2 h)}{60 h^4}-\frac{2 f(x+3 h)}{5 h^4}+\frac{7 f(x+4 h)}{240 h^4}=f^{(4)}(x)+\frac{41 h^6 f^{(10)}(x)}{7560}+O\left(h^7\right) \]

五阶导

七点（二阶）

\[ -\frac{f(x-3 h)}{2 h^5}+\frac{2 f(x-2 h)}{h^5}-\frac{5 f(x-h)}{2 h^5}+\frac{5 f(x+h)}{2 h^5}-\frac{2 f(x+2 h)}{h^5}+\frac{f(x+3 h)}{2 h^5}=f^{(5)}(x)+\frac{1}{3} h^2 f^{(7)}(x)+O\left(h^3\right) \]

九点（四阶）

\[ \frac{f(x-4 h)}{6 h^5}-\frac{3 f(x-3 h)}{2 h^5}+\frac{13 f(x-2 h)}{3 h^5}-\frac{29 f(x-h)}{6 h^5}+\frac{29 f(x+h)}{6 h^5}-\frac{13 f(x+2 h)}{3 h^5}+\frac{3 f(x+3 h)}{2 h^5}-\frac{f(x+4 h)}{6 h^5}=f^{(5)}(x)-\frac{13}{144} h^4 f^{(9)}(x)+O\left(h^5\right) \]

十一点（六阶）

\[ -\frac{13 f(x-5 h)}{288 h^5}+\frac{19 f(x-4 h)}{36 h^5}-\frac{87 f(x-3 h)}{32 h^5}+\frac{13 f(x-2 h)}{2 h^5}-\frac{323 f(x-h)}{48 h^5}+\frac{323 f(x+h)}{48 h^5}-\frac{13 f(x+2 h)}{2 h^5}+\frac{87 f(x+3 h)}{32 h^5}-\frac{19 f(x+4 h)}{36 h^5}+\frac{13 f(x+5 h)}{288 h^5}=f^{(5)}(x)+\frac{139 h^6 f^{(11)}(x)}{6048}+O\left(h^7\right) \]

六阶导

七点（二阶）

\[ -\frac{20 f(x)}{h^6}+\frac{f(x-3 h)}{h^6}-\frac{6 f(x-2 h)}{h^6}+\frac{15 f(x-h)}{h^6}+\frac{15 f(x+h)}{h^6}-\frac{6 f(x+2 h)}{h^6}+\frac{f(x+3 h)}{h^6}=f^{(6)}(x)+\frac{1}{4} h^2 f^{(8)}(x)+O\left(h^3\right) \]

九点（四阶）

\[ -\frac{75 f(x)}{2 h^6}-\frac{f(x-4 h)}{4 h^6}+\frac{3 f(x-3 h)}{h^6}-\frac{13 f(x-2 h)}{h^6}+\frac{29 f(x-h)}{h^6}+\frac{29 f(x+h)}{h^6}-\frac{13 f(x+2 h)}{h^6}+\frac{3 f(x+3 h)}{h^6}-\frac{f(x+4 h)}{4 h^6}=f^{(6)}(x)-\frac{13}{240} h^4 f^{(10)}(x)+O\left(h^5\right) \]

十一点（六阶）

\[ -\frac{1023 f(x)}{20 h^6}+\frac{13 f(x-5 h)}{240 h^6}-\frac{19 f(x-4 h)}{24 h^6}+\frac{87 f(x-3 h)}{16 h^6}-\frac{39 f(x-2 h)}{2 h^6}+\frac{323 f(x-h)}{8 h^6}+\frac{323 f(x+h)}{8 h^6}-\frac{39 f(x+2 h)}{2 h^6}+\frac{87 f(x+3 h)}{16 h^6}-\frac{19 f(x+4 h)}{24 h^6}+\frac{13 f(x+5 h)}{240 h^6}=f^{(6)}(x)+\frac{139 h^6 f^{(12)}(x)}{12096}+O\left(h^7\right) \]