Mathematica 中微分方程离散化形式非常简单，我们先列出 Mathematica 中的差分形式。网页 数值直线方法，或者 Mathematica 的文档页 tutorial/NDSolveMethodOfLines 中可以找到。

一阶差分

二阶差分

如果列表中没有，但是又需要使用 n 阶 m 点差分形式的时候怎么办呢？这就需要使用我们上次在《微分方程离散化》中讲到的待定系数法，Mathematica 中实现非常简单。

待定系数法生成差分形式

下面以 3 阶导的 9 点差分形式为例：

In[]:= 
1
n = 3

表示 3 阶导。

In[]:= 
1
tab = {-5, -4, -3, -2, 0, 1, 2, 3, 4}

表示使用了 附近的 9 个点：

系数分别为

注意这里作为演示特意去掉了其中一个点 。

In[]:= 
1
2
m = Length[tab];
partial = Table[0, n]~Join~{1}~Join~Table[0, m - n - 1]
Out[]:=

m 表示使用的是 9 点差分形式。第二行表示求待定系数的时候，除了第三阶导数的系数为 1，其它均为 0，需要注意的是第一个表示 0 阶导，也就是 自身。

In[]:= 
1
sum = Sum[a[i] Series[f[x + i h], {h, 0, m + 1}], {i, tab}]
Out[]:=

方程右边做泰勒展开，展开阶数通常比点数多 2，最后求误差项的时候需要用到。点数加 1 那一项有可能为 0，所以就需要加 2。

In[]:= 
1
2
3
sol = Solve[partial == Table[Coefficient[sum, \!\(\*SuperscriptBox[\(f\), 
TagBox[RowBox[{"(", "j", ")"}], Derivative],
MultilineFunction->None]\)[x]], {j, 0, m - 1}], Table[a[i], {i, tab}]]
Out[]:=

求解出所有系数。输出差分形式:

In[]:= 
1
Sum[a[i] f[x + i h], {i, tab}] /. sol[[1]] // Simplify
Out[]:=

下面输出误差项：

In[]:= 
1
2
3
Coefficient[sum, \!\(\*SuperscriptBox[\(f\), 
TagBox[RowBox[{"(", "m", ")"}], Derivative],
MultilineFunction->None]\)[x]] /. sol[[1]]
Out[]:=

如果这一阶误差为 0，那么可以输出下一阶的误差：

In[]:= 
1
2
3
Coefficient[sum, \!\(\*SuperscriptBox[\(f\), 
TagBox[RowBox[{"(", RowBox[{"m", "+", "1"}], ")"}], Derivative],
MultilineFunction->None]\)[x]] /. sol[[1]]
Out[]:=

如此绝妙的算法不收藏、点赞、转发一下吗？😄

NDSolve`FiniteDifferenceDerivative[]

Mathematica 中的 ‘NDSolve`FiniteDifferenceDerivative[]’ 函数可以用来生成差分格式，因此在数值求解微分方程的时候非常有用。

In[]:= 
1
2
3
xGrid = h Range[0, 4];
NDSolve`FiniteDifferenceDerivative[1, xGrid, Outer[f, xGrid], 
 "DifferenceOrder" -> 2]
Out[]:=

DifferenceOrder 表示采用几点进行差分，中间段采用该设置，两端采用更多点的差分形式。

In[]:= 
1
2
3
4
xGrid = h Range[0, 3];
yGrid = h Range[0, 4];
NDSolve`FiniteDifferenceDerivative[{1, 2}, {xGrid, yGrid}, 
  Outer[f, xGrid, yGrid], "DifferenceOrder" -> 2] // Simplify

表示对 x 方向一阶导，y 方向二阶导。x 方向和 y 方向的间隔 h 可以不同。