上次在《不动点迭代法数值求解泊松方程》和《不动点迭代法数值求解热传导方程》中讲到不动点迭代法效率特别低，这次讲在 Mathematica 中使用隐式法数值求解微分方程。该方法是一种通用数值算法，非线性方程组、常微分方程组、偏微分方程组、甚至刚性问题，无论是初值问题（initial value problem）还是边值问题（boundary value problem），都可以使用隐式法求解。如此精妙的算法我都舍不得分享😂

隐式法数值求解一维泊松方程

还是考虑普遍的形式：

并考虑边值条件：

将 x 的范围从 0 到 1 分割成 100 份：

In[]:= 
1
2
3
h = 0.01;
n = Floor[1/h]
xGrid = Range[0, 1, h];

边值条件，以及使用前面讲的离散化方法写出差分方程组：

In[]:= 
1
2
3
4
5
u[0] = 0; u[n] = 0;

equations = 
  Table[(u[i + 1] - 2 u[i] + u[i - 1])/h^2 + (u[i + 1] - u[i - 1])/(
     2 h) + u[i] + (i h) == 0, {i, 1, n - 1}];

该方程组为非线性方程组，下面讲两种求解方法。第一种利用 NSolve[] 函数：

In[]:= 
1
NSolve[equations, Map[u, Range[1, n - 1]]]

上面方程组不太难，所以可以使用 NSolve[] 函数求解。但是对于高度非线性的方程组来说，NSolve[] 无法给出解，所以下面讲第二种方法，利用 FindRoot[] 函数求解。FindRoot[] 函数需要猜测一个初始值，下次讲到牛顿法（Newton’s method）就知道为什么需要初始值。

In[]:= 
1
2
3
try = Table[{u[i], 0.5}, {i, 1, n - 1}];

sol = FindRoot[equations, try]
Out[]:=

下面画出了解析解和数值解的对比:

In[]:= 
1
2
Show[{ListPlot[Transpose[{xGrid, Table[u[i], {i, 0, n}] /. sol}]],
  Plot[uSol[x], {x, 0, 1}, PlotStyle -> Red]}]
Out[]:=

隐式法数值求解二维泊松方程

考虑普遍的形式：

并考虑狄利克雷边界条件：

离散求解区域：

In[]:= 
1
2
3
4
5
6
7
8
dx = 0.02;
dy = 0.04;

xGrid = Range[-2, 2, dx];
yGrid = Range[-2, 2, dy];

nx = Length[xGrid];
ny = Length[yGrid];

边值条件，以及离散化形成差分方程组：

In[]:= 
1
2
3
4
5
6
equations = Flatten[
  Table[If[xGrid[[i]]^2 + yGrid[[j]]^2 <= 1, u[i, j] == 1, 
    If[xGrid[[i]]^2 + yGrid[[j]]^2 >= 4, u[i, j] == 0, 
      (u[i + 1, j] - 2 u[i, j] + u[i - 1, j])/dx^2 + (
       u[i, j + 1] - 2 u[i, j] + u[i, j - 1])/dy^2 == u[i, j]]], 
       {i, 1, nx}, {j, 1, ny}], 1]

利用 FindRoot[] 函数求解：

In[]:= 
1
2
3
try = Flatten[Table[{u[i, j], 0.5}, {i, 1, nx}, {j, 1, ny}], 1];

sol = FindRoot[equations, try]
Out[]:=

下面画出数值解：

In[]:= 
1
2
3
4
uuSol = Flatten[Table[{xGrid[[i]], yGrid[[j]], u[i, j] /. sol}, {i, nx}, {j, ny}], 1];

ListDensityPlot[uuSol, PlotRange -> All, 
 RegionFunction -> Function[{x, y, z}, 1 < Norm[{x, y}] < 2]]
Out[]:=

牛顿法

FindRoot[] 函数内置了很多数值求解方程组的方法，可以使用 Method-> 指定。我们前面已经讲了二分法和不动点迭代法，下次介绍牛顿法（Newton’s method）。

请大家多多关注、点赞、收藏、转发，这一章绝对干货满满，以后会分享更多实用的算法。

万分感谢！🙏