上次在《牛顿法数值求解泊松方程》中讲了使用牛顿法（Newton’s method）数值求解多变量方程，并以一维泊松方程为例进行隐式求解。本次讲解使用牛顿法数值求解多维方程。迭代形式为上次所讲：

还是考虑方程：

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

离散求解区域：

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
f = {};
var = {};
index = {};

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

In[]:= 
1
2
3
4
5
6
Do[If[xGrid[[i]]^2 + yGrid[[j]]^2 <= 1, u[i, j] = 1, 
  If[xGrid[[i]]^2 + yGrid[[j]]^2 >= 4, u[i, j] = 0, 
  (AppendTo[f, (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]];
    AppendTo[var, u[i, j]]; AppendTo[index, {i, j}])
   ]], {i, 1, nx}, {j, 1, ny}]

雅可比行列式为：

In[]:= 
1
jacobian = D[f, {var}];

初始试探解为：

In[]:= 
1
2
3
4
uu = Table[0.5, {i, 1, Length[var]}];
Do[u[index[[i, 1]], index[[i, 2]]] = uu[[i]], {i, 1, Length[var]}];
error = Abs[Total[f]]
nn = 0;

下面使用牛顿法（Newton’s method）求解该问题：

In[]:= 
1
2
3
4
5
6
7
While[error > 10^-10 && nn < 100,
 uu = uu - LinearSolve[jacobian, f];
 
 Do[u[index[[i, 1]], index[[i, 2]]] = uu[[i]], {i, 1, Length[var]}];
 error = Abs[Total[f]];
 nn++;
 ]

迭代次数：

In[]:= 
1
nn
Out[]:=

误差：

In[]:= 
1
error
Out[]:=

下面画出数值解：

In[]:= 
1
2
uuSol = Flatten[Table[{xGrid[[i]], yGrid[[j]], u[i, j]}, {i, nx}, {j, ny}], 1];
ListDensityPlot[uuSol, RegionFunction -> Function[{x, y, z}, 1 < Norm[{x, y}] < 2]]
Out[]:=

请大家多多关注、点赞、收藏、转发，让更多有需要的人可以看到，以后会分享更多实用的算法。

万分感谢！🙏