数值分析·学习 | 平方根法和追赶法matlab实现

一、前言：

个人学习内容分享

二、算法描述：

平方根法：

如果A为n阶对称正定矩阵，则存在一个实的非奇异下三角矩阵L,使 $A=LL^{T}$ 编辑,当限定L的对角元素为正时，这种分解是唯一的。

下面我们用直接分解方法来确定计算L元素的递推公式，因为

$A=\begin{pmatrix} l_{11} & 0 & \cdots &0 \\ l_{21}&l_{22} & \cdots &0 \\ \vdots &\vdots &\ddots &\vdots \\ l_{n1}&l_{n2} & \cdots &l_{nn} \end{pmatrix}\begin{pmatrix} l_{11} & l_{21} & \cdots & l_{n1} \\ 0&l_{22} & \cdots & l_{n2} \\ \vdots &\vdots &\ddots &\vdots \\ 0&0 & \cdots &l_{nn} \end{pmatrix}$ 编辑

其中 $l_{ii}>0(i=1,2,\cdots ,n)$ 编辑.由矩阵乘法及 $l_{jk}=0$ 编辑(当j<k时)，得

$A_{ij}=\sum_{k=1}^{n}l_{ik}l_{jk}=\sum_{k=1}^{j-1}l_{jk}l_{jk}+l_{ij}l_{ij}$ 编辑

于是得到对于求解Ax=b的平方根法计算公式：

对于j=1,2,...,n

$(1) \; l_{ij}=(a_{ij}-\sum_{k=1}^{j-1}l_{ij}^{2})^{1/2};$ 编辑

$(2)\; l_{ij}=\frac{a_{ij}-\sum_{k=1}^{j-1}l_{ik}l_{jk}}{l_{jj}};i=j+1,\cdots ,n.$ 编辑

求解Ax=b，即求解两个三角形方程组：

①Ly=b,求y; ②LTx=y,求x.

$(3)\; y_{i}=\frac{b_{i}-\sum_{k=1}^{i-1}l_{ik}y_{k}}{l_{ii}};i=1,2,\cdots ,n.$ 编辑

$(4)\; x_{i}=\frac{b_{i}-\sum_{k=i+1}^{n}l_{ki}x_{k}}{l_{ii}};i=n,n-1,\cdots ,1.$ 编辑

其中： $a_{jj}=\sum_{k=1}^{j}l_{jk}^2;i=1,2,\cdots ,n.$ 编辑

$l_{jk}^2\leq a_{jj}\leq \underset{1\leq j \leq n}{max}\left \{a_{jj} \right \}$ 编辑

所以 $max_{jk}\left \{ a_{jk}^2 \right \}\leq\underset{1\leq j \leq n}{max}\left \{ a_{jj} \right \}$ 编辑

改进的平方根法：

为了避免开方，我们下面用定理9的分解式A= $LDL^{T}$ 编辑

即 $A=\begin{pmatrix} 1 &0 & \cdots &0 \\ l_{21}&1 & \cdots & 0\\ \vdots & \vdots & \ddots &\vdots \\ l_{n1}& l_{n2} & \cdots &1 \end{pmatrix}\begin{pmatrix} d_1 &0 & \cdots &0 \\ 0&d_2 & \cdots & 0\\ \vdots & \vdots & \ddots &\vdots \\ 0& 0 & \cdots &d_n \end{pmatrix}\begin{pmatrix} 1 &l_{21} & \cdots &l_{n1} \\ 0&1 & \cdots & l_{n2}\\ \vdots & \vdots & \ddots &\vdots \\ 0& 0 & \cdots &1 \end{pmatrix}$ 编辑

由矩阵乘法，并注意 $l_{jj}=1,l_{jk}=0(j<k)$ 编辑,得

$a_{ij}=\sum_{k-1}^{n}(LD)_{ik}(LT)_{kj}=\sum_{k=1}^{n}l_{ik}d_{k}l_{jk}=\sum_{k=1}^{j-1}l_{ik}d_{k}l_{jk}+l_{ij}d_{j}l_{jj}$ 编辑

于是得到计算L得元及D得对角元素公式：

对于i=1,2,...,n.

$(1)\; l_{ij}=\frac{a_{ij}-\sum_{k=1}^{j-1}l_{ik}d_{k}l_{jk}}{d_{j}},j=1,2,\cdots ,i-1;$ 编辑

$(2)\; d_{i}=a_{ii}-\sum_{k=1}^{j-1}l_{ik}^{2}d_{k}.$ 编辑

为了避免重复计算，我们引进tij=lijdj,由上式得到按行计算L,T元素得公式：d1=a11

对于i=2,3,...,n.

$(1)\; t_{ij}=a_{ij}-\sum_{k=1}^{j-1}t_{ik}l_{jk},j=1,2,\cdots ,i-1;$ 编辑

$(2)\; l_{ij}=\frac{t_{ij}}{d_{j}},j=1,2,\cdots ,i-1;$ 编辑

$(3)\; d_{i}=a_{ii}-\sum_{k=1}^{j-1}t_{ik}l_{jk}.$ 编辑

计算出T=LD的第i行元素tij(j=1,2,...,i-1)后，存放在A的第i行相应位置，然后再计算L的第i行元素，存放在A的第i行.D的对角线元素存放在A的相应位置。例如

$A=\begin{pmatrix} a_{11} & 0 &0 &0 \\ a_{21} & a_{22} &0 &0 \\ a_{31} & a_{32} & a_{33} &0 \\ a_{41} & a_{42} & a_{43} &a_{44} \end{pmatrix}\rightarrow \begin{pmatrix} d_{1} & 0 &0 &0 \\ l_{21} & d_{2} &0 &0 \\ l_{31} & a_{32} & d_{3} &0 \\ l_{41} & a_{42} & a_{43} &d_{4} \end{pmatrix}$ 编辑

对称正定矩阵A按 $LDL^{T}$ 编辑分解和按 $LL^{T}$ 编辑分解计算量差不多，但 $LDL^{T}$ 编辑分解不需要开方计算。

求解Ly=b,DLTx=y计算公式：

$(3)\begin{cases} & \text{ } y_{1}= b_{1}\\ & \text{ } y_{i}=b_{i}-\sum_{k=1}^{i-1} l_{ik}y_{k},i=2,3,\cdots ,n. \end{cases}$ 编辑

$(4)\begin{cases} & \text{ } x_{n}= \frac{y_{n}}{d_{n}}\\ & \text{ } x_{i}=\frac{y_{n}}{d_{n}}-\sum_{k=i+1}^{n}l_{ki}x_{k},i=n-1,\cdots ,2,1. \end{cases}$ 编辑

追赶法：

对于A为三对角矩阵时，可分解成以下形式

$A=\begin{pmatrix} b_{1} & c_{1} & \cdots & \cdots & 0 &0 \\ a_{2} & b_{2} & \cdots & \cdots & 0 &0 \\ \vdots & \vdots & \ddots & \ddots & \vdots &\vdots \\ \vdots & \vdots & \ddots & \ddots & \vdots &\vdots \\ 0 & 0 & \cdots & \cdots & b_{n-1} &c_{n-1} \\ 0 & 0 & \cdots & \cdots & a_{n} &b_{n} \end{pmatrix}= \begin{pmatrix} a_{1} & \cdots & \cdots &0 \\ r_{2} & \ddots & \cdots &0 \\ \vdots & \ddots & \ddots & \vdots \\ 0 & \cdots & r_{n} &a_{n} \end{pmatrix}\begin{pmatrix} 1 & B_{1} & \cdots &0 \\ 0 & 1 & \cdots &0 \\ \vdots & \vdots & \ddots &\vdots \\ 0 & 0 & \cdots &B_{n-1} \\ 0 & 0 & \cdots &1 \end{pmatrix}$ 编辑

即可得到

$\begin{cases} & \text{ } b_{1}=a_{1},c_{1}=a_{1}B_{1} \\ & \text{ } a_{1}=r_{1},b_{i}=r{i}B_{i-1}+a_{i},i=2,3,\cdots,n \\ & \text{ } c_{i}=a_{i}B_{i},i=2,3,\cdots,n-1. \end{cases}$ 编辑

从而求解Ax=f等价于解两个三角形方程组：

①Ly=f,求y,②Ux=y,求x.

进而追赶法公式为：

（1）计算{Bi}的递推公式

$B_{1}=c_{1}/b_{1};$ 编辑

$B_{i}=c_{i}/(b_{i}-a_{i}B_{i-1}),i=2,3,\cdots,n-1;$ 编辑

（2）解Ly=f

$y_{1}=f_{1}/b_{1},$ 编辑

$y_{i}=(f_{i}-a_{i}y_{i-1})/(b_{i}-a_{i}B_{i-1}),i=2,3,\cdots,n-1;$ 编辑

（3）解Ux=y

$x_{n}=y_{n},$ 编辑

$x_{i}=y_{i}-B_{i}x_{i+1},i=n-1,n-2,\cdots,2,1.$ 编辑

三、实现代码：

1、平方根法：

function X=squareRoot(A,b)
%%功能：平方根算法
n=length(A);
y=b;
X=y;
%变换矩阵
for k=1:n
    A(k,k)=sqrt(A(k,k));
    A(k+1:n,k)=A(k+1:n,k)/A(k,k);
    for j=k+1:n
        A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);
    end
end
%求解矩阵
for k=1:n
    y(k)=(b(k)- sum(A(k,1:k-1)*y(1:k-1)))/A(k,k);
end
A=A';
for k=n:-1:1
    X(k)=(y(k)-sum(A(k,k+1:n)*X(k+1:n)))/A(k,k);
end

2、改进的平方根法:

function X=ImprovedSquareRoot(A,b)
%%功能:改进的平方根公式
n=length(A);
v=b;
y=b;
X=b;
%变换矩阵
for j=1:n
    for k=1:j-1
        v(k)=A(j,k)*A(k,k);
    end
    A(j,j)=A(j,j)-sum(A(j,1:j-1)*v(1:j-1));
    for i=j+1:n
        A(i,j)=(A(i,j)-sum(A(i,k)*v(1:j-1)))/A(j,j);
    end
end
%求解矩阵
D=diag(diag(A));
A=A-D+diag(diag(ones(3)));
for k=1:n
    y(k)=(b(k)-sum(A(k,1:k-1)*y(1:k-1)))/A(k,k);
end
d=diag(D);
X(n)=y(n)/d(n);
for k=n-1:-1:1
    X(k)=y(k)/d(k)-sum(A(k+1:n,k).*X(k+1:n));
end
end

3、追赶法:

function X=CatchUpMethod(a,b,c,f)
%%功能：追赶法
n=length(b);
B=f;
B(1)=c(1)/b(1);
for k=2:n-1
    B(k)=c(k)/(b(k)-a(k)*B(k-1));
end
y(1)=f(1)/b(1);
for k=2:n
    y(k)=(f(k)-a(k-1)*y(k-1))/(b(k)-a(k-1)*B(k-1));
end
X(n)=y(n);
for k=n-1:-1:1
    X(k)=y(k)-B(k)*X(k+1);
end
end