上次在《牛顿法数值求解二维热传导方程》中使用牛顿法（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
68
69
70
71
72
73
74
75
76
(*Mathematica自身求解*)
\[CapitalOmega] = ImplicitRegion[1 <= x^2 + y^2 <= 4, {x, y}];

uSol = NDSolveValue[{D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] == u[x, y],
   DirichletCondition[u[x, y] == 1, x^2 + y^2 == 1],
   DirichletCondition[u[x, y] == 0, x^2 + y^2 == 4]},
  u, Element[{x, y}, \[CapitalOmega]]]

Plot3D[uSol[x, y], Element[{x, y}, \[CapitalOmega]]]

(*共轭梯度法*)
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]]

(*离散化*)
dx = 0.04;
dy = 0.05;

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

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

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


(*雅可比矩阵*)
jacobian = D[f, {var}];

(*初始试探解*)
uu = Table[0.5, {i, 1, Length[var]}];

(*牛顿法*)
error = 1; nn = 0;

While[error > 10^-8 && nn < 10000,
 
 Do[u[index[[i, 1]], index[[i, 2]]] = uu[[i]], {i, 1, Length[var]}];
 jaco = jacobian;
 jacoTrans = ConjugateTranspose[jaco];
 ff = f;
 error = Norm[ff];
 
 uu -= ConjugateGradient[jacoTrans . jaco, jacoTrans . ff];
 
 nn++]

输出：

In[]:= 
1
2
nn
error

画图：

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

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

万分感谢！🙏