上次在《牛顿法数值求解泊松方程》中使用牛顿法（Newton’s method）数值求解泊松方程。思路是将微分方程离散化成非线性方程组，然后使用牛顿法将非线性方程组简化成求解线性方程组，之后又介绍了《共轭梯度法数值求解线性方程组》。有朋友尝试将两者结合，发现求解出了问题，本次就附上完整代码以证明之前的牛顿法和共轭梯度法都没有任何问题。

牛顿法+共轭梯度法数值求解泊松方程

In[]:= 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
(*Mathematica自身求解*)
uSol[x_] = NDSolveValue[
  {u''[x] + u'[x] + u[x]^2 + x == 0,
   u[0] == 0, u[1] == 0},
  u[x], {x, 0, 1}]

(*共轭梯度法*)
ConjugateGradient[A_, b_] :=
 Module[{x, r, rr, p, Ap, k, \[Alpha], \[Beta]},
  
  x[0] = b;
  p[0] = r[0] = b - A . x[0];
  k = 0;
  
  While[k < 10000,
   rr = ConjugateTranspose[r[k]] . r[k];
   Ap = A . p[k];
   
   \[Alpha][k] = rr/ConjugateTranspose[p[k]] . Ap;
   x[k + 1] = x[k] + \[Alpha][k]  p[k];
   r[k + 1] = r[k] - \[Alpha][k]  Ap;
   
   (*Print[Norm[r[k+1]],",",x[k+1]];*)
   If[Norm[r[k + 1]] < 10^-10, Break[]];
   
   \[Beta][k] = ConjugateTranspose[r[k + 1]] . r[k + 1]/rr;
   p[k + 1] = r[k + 1] + \[Beta][k]  p[k];
   
   k++];
  x[k + 1]]

(*离散化*)
h = 0.01;
n = Floor[1/h]
xGrid = Range[0, 1, h];

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]^2 + (i  h),
   {i, 1, n - 1}];


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

(*雅可比行列式*)
jacobian = D[f, {var}];

(*初始试探解*)
uu = Table[0., {i, 0, n}];
uu[[1]] = 0; uu[[-1]] = 0;

(*牛顿法*)
error = 1; nn = 0;
While[error > 10^-8 && nn < 10000,
 
 Do[u[i] = uu[[i + 1]], {i, 1, n - 1}];
 jaco = jacobian;
 jacoTrans = ConjugateTranspose[jacobian];
 ff = f;
 error = Norm[ff];
 
 uu[[2 ;; -2]] -= 
  ConjugateGradient[jacoTrans . jaco, jacoTrans . ff];
 
 nn++]

代码中很多看似无用的设计，其实是为了加快运算速度。代码也很容易用其它语言实现，以后如果有机会了会发布 Python 或者 C++ 代码。

输出：

In[]:= 
1
2
3
nn
error
uu

画图：

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

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

万分感谢！🙏