上次在《牛顿法数值求解方程》中讲了使用牛顿法（Newton’s method）数值求解单变量方程。这次使用牛顿法数值求解多变量方程，并以泊松方程为例进行隐式求解，参考《隐式法数值求解泊松方程》。该方法极其实用，如此精妙的算法我都舍不得分享😂

使用牛顿法数值求解多变量方程无非是将单变量形式：

变成：

其中 为未知变量组成的向量， 为方程不为零的那边组成的向量，

为雅可比矩阵。

为求 ，需要求解线性方程组 ，其中 为未知向量。解此线性方程组可以使用共轭梯度法（Conjugate gradient method），以后有机会会介绍此方法。Mathematica 中求线性方程组的方法为 LinearSolve[]。

牛顿法数值求解一维泊松方程

还是考虑普遍的形式：

并考虑边值条件：

将 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
u[0] = 0; u[n] = 0;

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

变量为：

In[]:= 
1
var = Table[u[i], {i, 1, n - 1}];

雅可比行列式为：

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

初始试探解为：

In[]:= 
1
2
uu = Table[0.5, {i, 0, n}];
uu[[1]] = 0; uu[[-1]] = 0;

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

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

迭代次数：

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

误差：

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

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

In[]:= 
1
2
Show[{ListPlot[Transpose[{xGrid, uu}]],
  Plot[uSol[x], {x, 0, 1}, PlotStyle -> Red]}]
Out[]:=

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

万分感谢！🙏