Giáo trình Matlab căn bản - Chương 3: Hệ phương trình đại số tuyến tính

 end 
end 
i 
function x = pythag(y,z) 
%tinh sqrt( y^2 + z^2 ). 
rmax = max(abs([y;z])); 
if (rmax == 0) 
  x = 0; 
else 
  x = rmax*sqrt((y/rmax)^2 + (z/rmax)^2); 
end 
Để giải hệ phương trình ta dùng chương trình ctlsqr.m: 
clear all, clc 
maxiter = 50; 
A = [ 1 3 4; 2 5 7; 3 1 2]; 
b = [8  14  6]ʹ; 
x = lsqr(A, b, maxiter) 
§26. PHƯƠNG PHÁP SYMMLQ 
Liên quan đến phương pháp MINRES và CG là thuật toán SYMMLQ do 
Paige và Saunders đưa ra. Ta xét hệ phương trình [A][X] = [B] với [A] là ma 
trận đối xứng nhưng không cần xác định dương. Ta chọn nghiệm ban đầu là 
β1[v1] = {B],  [ ]1 2Bβ = . Tại  lần  lặp  thứ k của phương pháp CG  ta có được xk 
sao cho [rk] = [B] ‐ [A][Xk] trực giao. Do [Vk] là cơ sở trực giao nên ta có thể đặt 
[Xk] = [Vk][yk] và có: 
[rk] = [B] ‐ [A][Vk][yk] = β1[v1] ‐ [Vk][Tk][yk] ‐  [ ] [ ] [ ]Tk 1 k k 1( e y ) v+ +β   (1) 
Do [ ] [ ]Tk kV r 0=  nên nhân (1) với [ ]TkV và dùng [ ] [ ]Tk k 1V v 0+ =  và [ ]Tk 1 1V v e= ta 
có: 
  [ ] [ ] [ ][ ]Tk k 1 1 k k0 V r e T y= = β −               (2) 
để giải hệ (2), Paige và Saunders đề nghị thực hiện phép phân tích LQ: 
  [ ] [ ][ ] [ ] [ ] [ ]T Tk k k k kT L Q Q Q E= =  
 201
với  [ ]kL   là  ma  trận  tam  giác  và  [ ]kQ   là  ma  trận  trực  giao.  Thuật  toán 
SYMMLQ gồm các bước sau: 
  ‐ Cho x0 
  ‐ Tính x = xo, r = b ‐ Ax,  rρ = , v r= ρ  
  ‐ β = 0,  0β =% , c = ‐1, s = 0, k = ρ 
  ‐ vold = 0, w = v, g = 0, g 0=%  
  ‐ Lặp khi k < tol: 
    •  oldv Av v= − β%  
•  v * vα = % , v v v= − α% %  
•  vβ = % ,  oldv v= , v v /= β%  
•  1l s c= α − β% ,  2l s= β  
•  s cα = − β − α%% ,  cβ = β%  
•  2 20l = α + β%  
•  0 0c l ,s l= α = β%  
•  1 2 0ˆ ˆg g l g, g l g, g g l= − = − =% %  
• x x (gc)w+(gs)v= +  
• w = sw ‐ cv  
•  2 2ˆk g g= +%  
Ta xây dựng hàm symmlq() để thực hiện thuật toán này: 
function x = symmlq(A, b, x, maxiter, tol) 
%Ham thuc hien thua toan SYMMLQ voi A la ma tran doi xung 
[m,n] = size(A); 
n2b = norm(b);                       
xmin = x;                          
imin = 0;                          
tolb = tol * n2b;               
r = b ‐ A * x;                 
normr = norm(r);               
normrmin = normr;                   
v = r; 
vold = r; 
u = vold; 
v = u; 
 202
beta1 = voldʹ * v; 
beta1 = sqrt(beta1); 
vv = v / beta1; 
wbar = vv; 
v = A * vv; 
alpha = vvʹ * v; 
v = v ‐ (alpha/beta1) * vold; 
numer = vvʹ * v; 
denom = vvʹ * vv; 
v = v ‐ (numer/denom) * vv; 
volder = vold; 
vold = v; 
u = vold; 
v = u; 
betaold = beta1; 
beta = voldʹ * v; 
beta = sqrt(beta); 
gammabar = alpha; 
deltabar = beta; 
gamma = sqrt(gammabar^2 + beta^2); 
cs = gammabar / gamma; 
sn = beta / gamma; 
zeta = beta1 / gamma; 
epsilonzeta = 0; 
normrcgcs = abs(beta1 * sn); 
if (cs == 0) 
    normrcg = Inf; 
else 
    normrcg = normrcgcs / abs(cs); 
end 
stag = 0;                           
for i = 1 : maxiter    
   vv = v / beta; 
   w = cs * wbar + sn * vv; 
   stagtest = zeros(n, 1); 
   ind = (x ~= 0); 
 203
   stagtest(ind) = w(ind) ./ x(ind); 
   stagtest(~ind & (w ~= 0)) = Inf; 
   if (zeta == 0) | (abs(zeta)*norm(stagtest, inf) < eps) 
       stag = stag + 1; 
   else 
       stag = 0; 
   end 
   x = x + zeta * w; 
   wbar = sn * wbar ‐ cs * vv;     
   v = A * vv;   
   v = v ‐ (beta / betaold) * volder; 
   alpha = vvʹ * v; 
   v = v ‐ (alpha / beta) * vold; 
   volder = vold; 
   vold = v;  
   u = vold;   
   v = u; 
   betaold = beta; 
   beta = voldʹ * v; 
   if (beta < 0) 
      break 
   end 
   beta = sqrt(beta); 
   delta = cs * deltabar + sn * alpha; 
   deltazeta = ‐ delta * zeta; 
   gammabar = sn * deltabar ‐ cs * alpha; 
   epsilon = sn * beta; 
   deltabar = ‐ cs * beta; 
   gamma = sqrt(gammabar^2 + beta^2); 
   csold = cs; 
   snzeta = sn * zeta; 
   cs = gammabar / gamma; 
   sn = beta / gamma; 
   epszdelz = epsilonzeta + deltazeta; 
   epsilonzeta = ‐ epsilon * zeta; 
   zeta = epszdelz / gamma; 
 204
   mrcg = norm((csold*epszdelz/gammabar ‐ snzeta)*vold);  
   normr = sqrt(epszdelz^2 + epsilonzeta^2); 
   normrcgcs = normrcgcs * abs(sn); 
   if (cs == 0) 
       normrcg = Inf; 
   else 
       normrcg = normrcgcs / abs(cs); 
   end   
   if (normr <= tolb)      
       r = b ‐ A * x;    
       normr = norm(r); 
       if (normr <= tolb) 
           break 
       end 
   end 
   if (normrcg <= tolb)  
       xcg = x + (epszdelz/gammabar) * wbar; 
       r = b ‐ A * xcg;                       
       normrcg = norm(r); 
       if (normrcg <= tolb) 
         x = xcg; 
         break 
      end 
   end 
   if (stag >= 2)                 
      break 
   end 
   if (normr < normrmin)          
      normrmin = normr; 
      xmin = x; 
      imin = i; 
   end 
end                          
r = b ‐ A * x;               
normr = norm(r);                  
i 
 205
Để  giải  hệ  phương  trình  bằng  thuật  toán  SYMMLQ  ta  dùng  chương  trình 
ctsymmlq.m: 
clear all, clc 
A = [ 1  2  4  1;2  3  1  5;4  1  1  6;1  5  6  5]; 
b = [ 8  11  12  17]ʹ; 
maxiter = 50; 
x = [0  0  0]ʹ; 
tol = 1e‐12; 
x = symmlq(A, b, x, maxiter, tol) 
§27. PHƯƠNG PHÁP CHEBYSHEV 
  Tính hội  tụ  của phương pháp  lặp phụ  thuộc vào  tính  chất  của phổ  ‐ 
nghĩa  là  của  các giá  trị  riêng  ‐  của ma  trận  [A].  Để  cải  thiện  tính  chất này 
người ta thường biến đổi hệ phương trình tuyến tính bằng một phép biến đổi 
tuyến tính thích hợp. Quá trình này được gọi là thử trước(preconditioner). Ví 
dụ nếu ma trận [M] xấp xỉ ma trận hệ số [A] theo một cách nào đó, hệ được 
biến đổi 
  [M]‐1[A][X] = [M]‐1[B] 
sẽ có nghiệm như hệ phương trình [A][X] = [B] nhưng tính chất phổ của hệ số 
của ma trận [M]‐1[A] có thể thuận lợi hơn. 
Ta xét phương trình với [A] là ma trận đối xứng, xác định dương. Lúc 
đó phổ của ma trận [A] sẽ nằm trong đoạn [λmin, λmax] với λmin, λmax là các giá 
trị riêng lớn nhất và nhỏ nhất của [M]‐1[A]. Thuật toán tìm nghiệm là: 
‐ cho [X0], tính [R0] = [B] ‐ [A][X0] 
‐ chọn tham số α và c sao cho phổ của [A] nằm trên đoạn [d ‐ c, d + c] 
hay trong ellip có tiêu điểm d ± c không chứa gốc toạ độ 
và tính với n = 1, 2,..., n cho đến khi hội tụ: 
  [Z] = [M]‐1[R] 
  2
d
α =   [P] = [Z]      khi n = 1 
  [ ] [ ]2n 1n n n n n‐1
n
c P Z P
2 d
−α 1⎛ ⎞β = α = = [ ] + β⎜ ⎟ − β⎝ ⎠         khi n ≥ 2 
  [ ] [ ] [ ]n 1 n n nX X P+ = + α  
[Rn+1] = [Rn] ‐ αn[A][Pn] 
 206
Ta xây dựng hàm chebyiter() để thực hiện thuật toán trên: 
function  x = chebyiter ( A, x, b, M, maxiter, tol ) 
%  Cu phap  x = chebyiter ( A, x, b, M, maxiter, tol ) 
%  Dung pp lap Chebyshev de giai he pt  A*x = b. 
%  A(n, n) ma tran doi xung, xac dinh duong 
%  X(n), vec to nghiem ban dau 
%  B(n), ve phai 
%  M, ma tran preconditioner  
%  cho M bang mt don vi neu khong thu truoc 
%  X(n), nghiem 
if size(x, 1) == 1 
    x = xʹ; 
end 
r = b ‐ A * x; 
eigs = eig ( inv ( M ) * A ); 
eigmax = max ( eigs ); 
eigmin = min ( eigs ); 
c = ( eigmax ‐ eigmin ) / 2.0; 
d = ( eigmax + eigmin ) / 2.0; 
for i = 1 : maxiter 
    z =  M \ r; 
    if ( i == 1 ) 
        p = z; 
        alfa  = 2.0 / d;         
    else 
        beta = ( c * alfa / 2.0 )^2; 
        alfa = 1.0 / ( d ‐ beta ); 
        p = z + beta * p;         
    end 
    x  = x + alfa * p; 
    r = r ‐ alfa * A * p; 
    err  = norm ( r ); 
    if ( err <= tol  ) 
        break  
    end 
 207
end 
Ta dùng chương trình ctchebyiter.m để giải hệ phương trình:  
clear all, clc; 
    n = 10; 
    A = zeros ( n, n ); 
    for i = 1 : n 
      A(i, i) = 3.0; 
    end 
for i = 1 : n‐1 
         A(i, i + 1) = ‐1; 
    end 
    for i = 1 : n‐1 
         A(i + 1, i) = ‐1; 
    end 
    x = [1:n ]ʹ; 
    b = A * x; 
    x = ones ( n, 1 ); 
    M = 2.0 * eye ( n ); 
    maxiter = 50; 
    tol = 1e‐6; 
    y = chebyiter ( A, x, b, M, maxiter, tol ); 
    fprintf(ʹNghiem cua he phuong trinh\nʹ);     
fprintf(ʹ    %f\nʹ, y) 
§28. PHƯƠNG PHÁP QR 
Ta phân tích ma trận hệ số [A] thành: 
  [A] = [Q][R] 
Do   [Q]T[Q] = [E]  
nên: 
  [A][X] = [Q][R][X] = [B] 
  [Q]T[A][X] = [Q]T[Q][R][X] = [R][X] = [Q]T[B] 
Do [R] là ma trận tam giác trên nên ta tìm nghiệm dễ dàng. Ta xây dựng hàm 
givens() để thực hiện phép quay Givens: 
 208
function [c, s, r] = givens(x, y); 
%  tinh c, s, r sao cho  [c  s] [x] = [r] 
%  [‐s c] [y] =  [0] 
%  voi c*c + s*s = 1; 
if (y == 0) 
    c = 1;  
    s = 0;  
    r = x; 
else  
    if (abs(x) >= abs(y)) 
        t = y/x;  
        r = sqrt(1 + t*t); 
        c = 1/r; 
        s = t*c; 
        r = x*r; 
    else  
        t = x/y;  
        r = sqrt(1 + t*t); 
        s = 1/r; 
        c = t*s; 
        r = y*r; 
    end 
end  
Tiếp  theo  ta  xây  dựng  hàm  qrgivens()  thực  hiện  việc  tìm  nghiệm  của  hệ 
phương trình bằng thuật toán phân tích QR nhờ phép quay Givens: 
function x = qrgivens(A, b); 
[m, n] = size(A); 
tau = zeros(n, 1); 
%R = [A(1:n+1, :) b(1:n+1)]; 
R = [A(1:n, :) b(1:n)]; 
for j = 2:n 
    for i = 1:j‐1      
        [c, s, r] = givens(R(i, i), R(j, i)); 
        R(i, i) = r;  
 209
        R(j, i) = 0; 
        t = c*R(i, i+1:n+1) + s*R(j, i+1:n+1); 
        R(j, i+1:n+1) = ‐s*R(i, i+1:n+1) + c*R(j, i+1:n+1); 
        R(i, i+1:n+1) = t; 
    end 
end 
for k = n+2:m, 
    a = [A(k, :) b(k)]; 
    for i = 1:n+1 
        [c, s, r] = givens(R(i, i),a(i)); 
        R(i,i) = r;  
        a(i) = 0; 
        t = c*R(i, i+1:n+1) + s*a(i+1:n+1); 
        a(i+1:n+1) = ‐s*R(i, i+1:n+1) + c*a(i+1:n+1); 
        R(i, i+1:n+1) = t; 
    end 
end 
x = R(1:n, n+1); 
   for j = n:‐1:2 
       x(j) = x(j)/R(j, j); 
       x(1:j‐1) = x(1:j‐1) ‐ R(1:j‐1, j)*x(j); 
   end  
   x(1) = x(1)/R(1, 1); 
Để giải hệ phương trình ta dùng chương trình ctqrgivens.m: 
clear all, clc 
A = [1 2 ‐1;2 1 1; 1 1 3]; 
b = [2 4 5]ʹ; 
x = qrgivens(A, b)