exp3_2.m

使用matlab软件编制的计算程序

% exp3_2.m --- 矩阵分解命令的学习

function study_lu_and_chol

% ---------- LU分解 ----------
% [简介] 设A是非奇异矩阵,则有如下列主元的LU分解
%             PA = LU
%        其中P是置换矩阵(又称排列矩阵),L是单位下三角矩阵且|l(ij)|<=1,U是上三角矩阵
%        矩阵A如果有了上述分解则求解线性方程组 Ax = b 就等价于分别求解两个三角方程组
%            Ly = Pb  和  Ux=y
%        而求解三角方程组是非常容易的.用这种方法求解方程组与列主元的Gauss消去法是等价的
%        但 LU 分解还有其它优点. 参见 P65 例9 和 P70 实验课题(一)

A = [10  -7  0
     -3   2  6
      5  -1  5];
[L U P] = lu(A)

% 如果再给b,则下面是解两个三角方程组的程序. 参见P49(3-15)式
b = [7 4 6]';
[n n]=size(A);
x = zeros(n,1);

% 前代求解:Ly = Pb(解用x储存)
b = P*b;
x(1) = b(1);
for k = 2:n
    x(k) = b(k)-L(k,1:k-1)*x(1:k-1);
end

% 回代求解:Ux = y (这里先y=x,解仍用x储存)
x(n) = x(n)/U(n,n);
for k = n-1:-1:1
    x(k) = ( x(k)-U(k,k+1:n)*x(k+1:n) ) / U(k,k);
end
x

% ---------- Cholesky分解 ----------
% [简介] 设A是(实对称)正定矩阵,则有如下唯一的分解
%             A = R'*R
%        其中 R 是上三角矩阵且主对角元全为正

% 例
A = [4    -1     1
    -1  4.25  2.75
     1  2.75   3.5];
R = chol(A)

% [注1]  有了Chol分解,求解 Ax = b 又可转化为分别求解两个三角方程组
%        这是容易的,请同学们自己完成(常称为平方根法)
% [注2]  P51 的 chol 分解 A = LDL'
%        避免了前面分解中的开方运算(开方运算比四则运算速度是慢的)



% ******* 你的实验  *******

% ★【实验一】
% 首先完成下面Gauss列选主元的消去法程序,单独存为 gauss.m 文件(注意一定要与函数名相同)
% 可参考P43 图3-2
% 然后找一个例子调用此程序验证是否正确(调用方法同Matlab内部函数调用完全一样)

% function x = gauss(A,b)
% [n,n] = size(A);
% x = zeros(n,1);
% 
% Aug = [A,b];                         % 增广矩阵
% 
% for k = 1:n-1
%     [piv,r] = max(abs(Aug(k:n,k)));  % 找列主元所在子矩阵的行r
%     r = r + k - 1;                   % 列主元所在大矩阵的行
%     
%     if r>k
%         ?                            % 对Aug实施行交换(一行命令就可以了,怎么写?)
%     end
%     
%     if Aug(k,k)==0, error('?'), end  % 程序遇到error会中断执行并显示其中的提示内容
%     
%     % 把增广矩阵消元成为上三角
%     for p = k+1:n
%         mult = Aug(p,k)/Aug(k,k);    % 消元乘子
%         Aug(p,k:n+1) = ?;
%     end
% end 
% 
% % 解上三角方程组
% A = Aug(:,1:n); b = Aug(:,n+1);
% x(n) = ?;
% for k = n-1:-1:1
%     x(k) = ?;
% end



% 【实验二】编下面程序
%  (1) 追赶法(P45) 
%  (2) Cholesky分解法(P51)

%  下面是追赶法的程序你可以参考( 参见 P45 ),注意向量 a 的下标与书上不同

% function d = tridiag(a,b,c,d)
% % d = tridiag(a,b,c,d) --- 求解三对角线性方程组的追赶法( 参见 P45 )
% % 方程组的形式
% %         |b1  c1            | |x1|     |d1|
% %         |a1  b2  c2        | |x2|     |d2|
% %         |   a2   b3  c3    | |x3|  =  |d3|
% %         |        a3  b4  c4| |x4|     |d4|
% %         |            a4  b5| |x5|     |d5|
% % 输入: a,b,c,d 是四个向量 
% % 输出: d 方程组的解

% n = length(d);
% k = 1;
% d(k) = d(k)/b(k);
% c(k) = c(k)/b(k);
% for k = 2:n-1
%     b(k) = b(k) - a(k-1)*c(k-1);
% 	  c(k) = c(k)/b(k);
%     d(k) = ( d(k) - a(k-1)*d(k-1) )/b(k);
% end
% k = n;
% b(k) = b(k)-a(k-1)*c(k-1);
% d(k) = ( d(k) - a(k-1)*d(k-1) )/b(k);
% 
% % d(n) = d(n);
% for k = n-1:-1:1 
% 	d(k) = d(k) - c(k)*d(k+1);
% end