matlab 实现共轭梯度法

2744557306

这段代码是关于共轭梯度法（Conjugate Gradient Method）的实现，用于解决线性代数系统 (Ax = b)。具体而言，它使用了预条件共轭梯度法（Preconditioned Conjugate Gradient, PCG）来求解具有对称正定系数矩阵 (A) 的线性方程组。
以下是代码的一些关键部分的解释：
: {! V# U  s  o# e5 c+ C8 {
1.(A) 和 (b)：给定的线性系统的系数矩阵和右侧向量。
# @" S- `& ^6 S7 ^5 C8 F2.(w)：一个权重参数，用于调整共轭梯度法的收敛性。
/ B* ?  V7 Y) ~# E3.(D)：(A) 的对角矩阵。
4.(CL) 和 (CLZ)：分别是 (A) 的严格下三角和严格上三角。
, ]+ }2 c& J% H- G3 o% k! i4 U) V5.(L)：预条件矩阵，通过 (L = (D - w \cdot CL) \cdot D^{1/2} / \sqrt{w \cdot (2 - w)}) 计算得到。
. Y+ m  h; N  @6.(M)： (M = L \cdot L^T)，用于预条件化。
7.(C)： (C = D^{-1} \cdot CL)。
8.(u)、(v)：初始的近似解和共轭梯度法中的辅助向量。
9.(rw)： (rw = g - B \cdot v)，其中 (B = L^{-1} \cdot A \cdot (L^{-1})^T)。
10.接下来是 PCG 的主要迭代过程，其中计算了共轭梯度法的一系列参数，如 (af)、(r1)、(zw)、(z1)、(bt)、(p)、(q) 等。
11.最终，通过迭代过程得到近似解 (u)。

1. A=[5,-4,1,0;-4,6,-4,1;1,-4,6,-4;0,1,-4,5];
   / j! t8 M3 R- H. P7 s) L, y+ R3 P! o
2. b=[2,-1,-1,2]';% f/ [! B+ W, u  ~; k! N/ O
   
3. n=length(b);
   9 l1 }: f. r6 Y2 I
4. w=10;
   * S/ `7 l7 M+ N( [  M
5. D=diag(diag(A));' @/ b+ }2 r& y' X5 J# q7 X2 y3 c) T
   
6. CL=-triu(A,1);) y$ C+ X$ n0 R& H
   
7. CLZ=CL';
   . N$ l& v& p. n8 p
8. L=((D-w*CL)*D.^(1/2))/sqrt(w*(2-w));, N0 S  o* P\" b# d
   
9. M=L*L';7 A1 e/ H+ l5 ?/ h
   
10. C=inv(D)*CL;
    , x. B+ C4 Q2 B; l
11. u=[2,3,4,5]';
    ' C# [/ |7 l. }7 c( h1 N6 a
12. g=inv(L)*b;0 o3 u! Z' o& Q5 ], R, Z
    
13. B=inv(L)*A*inv(L)';
    \" r( r6 v; y8 ]3 E! O7 ?
14. v=L'*u;( t+ [4 u/ B& Z* d1 j$ Y) _
    
15. rw=g-B*v;
    : d' B+ U% `/ m9 D& L$ B, r4 A
16. r=L*rw;
    3 x0 t: d# P$ S! b  I6 U7 C
17. p=inv(M)*r;4 T/ J* s6 K+ A
    
18. z=p;5 n1 t) j0 `7 t; k+ Y
    
19. q=A*p;- W: s$ u3 _: Z9 X) u9 W
    
20. for i=1:50  u* n* R* w6 ^# P  ]. b
    
21.     af=r'*z/(p'*q);, s$ v; ?. E/ ]0 R: O
    
22.     u=u+af*p
    ! C\" b\" F9 Q% K\" L) n4 U
23.     r1=r-af*q;/ U' S3 }- ?/ l& ~
    
24.     zw=(eye(n)-w*C)*D.^(1/2)\(w*(w-2)*r1);' F\" Y8 m5 \1 u; K6 F' q
    
25.     z1=D.^(1/2)*(eye(n)-w*C')\zw;( o% S8 h, e7 w1 g. U% e
    
26.     bt=r1'*z1/(r'*z);  g; o- H) m1 o! `5 z2 E
    
27.     p=z1+bt*p;
    . K  d5 w0 z& T
28.     q=A*p;
    6 T5 C$ ]! G4 U
29. end9 I# d# ]8 X* M. B7 K5 D9 g- E
    
30.   % Boundary condition.
    ; E7 W1 R9 i# g2 j- v* k
31. % zw(:,1) = 0;1 A* g1 w1 W8 I% v4 ~
    
32.   %zw(:,n) = 0;; k5 Y: w! g$ t' Z$ }3 l/ d
    
33.   %z(:,1) = 0;
    $ q; J& r8 q: v2 c# i
34.   %z(:,n) = 0;# l( c6 q' \' G# D  }* C; n6 J
    
35.   %for i=2:10
    1 L( Z: ?8 @) z8 m
36.    %   for j=2:10
    2 p# u% Q, o& ]' Z( ]% N
37.     %  zw(i,j)=w*(zw(i,j-1)+z(i-1,j))/4+w*(2-w)*

复制代码

。
" G, \: i, Q- F# |/ J3 q
+ X5 l7 X- A+ L' |: b7 G* `0 h