首页 资讯 :科学和工程计算中的线性方程组系数矩阵(图)

:科学和工程计算中的线性方程组系数矩阵(图)

更新时间:2022-09-25 8:08:34 分类:资讯 浏览:30

在科学和工程计算中平方根公式,经常需要求解形如 Ax=b 的线性方程组,其中 A 是 ntimes m 矩阵,称为系数矩阵,b 是 n 维列向量,称为 a右手向量,x为要求解的m维列向量平方根公式,称为解向量。

在科学和工程的实际计算中,经常会遇到系数矩阵A是一个对称正定矩阵。如果

A=begin{} a_{11} \ a_{21} & a_{22} && 对称\ a_{31} & a_{32} & a_{33} \ vdots & vdots & vdots & ddots \ a_{n1} & a_{n2} & dots & a_{nn} end{} \

对于正定矩阵,有如下三角矩阵

L = begin{} l_{11} \ l_{21} & l_{22} & & bm 0\ l_{31} & l_{32} & l_{33} \ vdots & vdots & vdots & ddots \ l_{n1} & l_{n2} & l_{n3} & dots & l_{nn} end{} \

令 A = L cdot L^{T} 成立。如果 L 的主要对角元素取正值,则这种分解是唯一的。

直接展开矩阵关系A = LL^{T},有

begin{align} a_{11} &= l_{11}^{2} \ a_{21} &= l_{21}l_{11},quad a_{22} = l_{21} ^{2}+l_{22}^{2}\a_{31} &= l_{31}l_{11},quad a_{32} = l_{31}l_{21}+l_{32} l_{22},quad a_{33}=l_{31}^{2}+l_{32}^{2}+l_{33}^{2}\dots end{align}\

据此,矩阵 L 的元素 l_{11} to l_{21} to l_{22} to l_{31} to l_{32} to dots 可以通过行,计算公式为

begin{cases} l_{ij} &= (a_{ij} - sum{k=1}^{j-1}l_{ik}l_{jk}) / l_{jj}, 四边形 & j = 1, 2, dots, i-1 \ l_{ii} &= (a_{ii} - sum{k=1}^{i-1}l_{ik}^2)^frac{1}{2}, quad & i = 1, 2, dots, n \ end{cases} \

基于矩阵分解 A = LL^{T} ,对称正定方程组 Ax = b 可以简化为两个三角方程 Ly = b 和 L^Tx = y 来求解。由 Ly = b 即

begin{cases} l_{11}y_{1} &= b_1 \ l_{21}y_{2} + l_{22}y_{2} &= b_2 \ dots dots dots l_{n1}y_{1} + l_{n2}y_{2} + dots + l_{nn}y_{n} &= b_n end{cases} \

y_1 to y_2 to dots to y_n 可以按顺序计算:

y_i = (b_i - sum{k=1}^{i-1}l_{ik}y_{k})/l_{ii}, quad i = 1, 2, dots,n

而 L^Tx = y 就是

begin{cases} begin{*}{2} l_{11}x_1 + l_{21}x_2 + dots + l_{n1}x_n &= y_1 \ l_{22}x_2 + dots + l_{n2}x_n &= y_2 \ dots dots dots \ l_{nn}x_n &= y_n \ end{*} end{cases} \

找到 x_n to x_{n-1} to dots to x_1 的可逆顺序:

x_i = (y_i - sum{k=i+1}^nl_{ki}x_{k})/l_{ii}, quad i = n, n-1, dots, 1 \

由于公式中包含矩阵分解时的平方根运算,所以该算法称为平方根法,也称为分解法。

根据上面的公式,可以通过编写程序求解方程:

subroutine cholesky_full(n, a, y)
    implicit none
    
    integer, intent(in) :: n
    real, intent(inout) :: a(n, n), y(n)
    
    integer :: i, j, k
    real :: temp
    
    ! 分解矩阵,生成下三角阵L
    ! 工程问题中的很多矩阵非常庞大,所以,计算过程中的数据应该直接存放在原始数组a中,
    ! 而不是新创建一个数组
    do i = 1, n
        ! 公式中,j的取值范围为1到j-1,此处换成1到j,可以将分解式统一起来,省去一次判断。
        ! 因为j=i时,j循环虽然会执行错误的操作、生成错误的a(i,j)结果,
        ! 但a(i,j)马上就会被最外层的i循环生成的正确数据替换
        do j = 1, i
            temp = a(i, j)
            do k = 1, j-1

1的平方加到100的平方公式_求x平方-2x 1=0的根_平方根公式

temp = temp - a(i, k) * a(j, k) end do a(i, j) = temp / a(j, j) ! a(j, i) = 0. ! 对角线上方d的元素赋0,可有可无 end do a(i, i) = sqrt(temp) end do ! 根据Ly=b求解出y do i = 1, n temp = y(i) do j = 1, i-1 temp = temp - a(i, j) * y(j) end do y(i) = temp / a(i, i) end do ! 求解出x do i = n, 1, -1 temp = y(i) / a(i, i) y(i) = temp ! 公式中k的范围为i+1到n,此处为1到i-1,因为下方a(i,k)的下标和公式中交换了顺序 do k = 1, i-1 y(k) = y(k) - temp * a(i, k) end do end do end subroutine

上面代码的分解部分和前面的公式基本一样,容易理解,只是引入了一个临时变量temp来存储数据。但是,如果我们交换j和k两层循环的位置,然后稍微调整循环计数器的取值范围,就可以直接完成分解操作,而无需使用临时变量。代码如下:

do i = 1, n
    do k = 1, i - 1
        a(i, k) = a(i, k) / a(k, k)
        do j = k + 1, i
            a(i, j) = a(i, j) - a(i, k) * a(j, k)
        end do
    end do
    a(i, i) = sqrt(a(i, i))
end do

参考资料:

王能超。高等学校教材,数值分析简明教程,(第2版)[M]. 2003.

吴建平、王正华、李晓梅。稀疏线性方程组的高效求解与并行计算[M].湖南科学技术出版社,2004.

版权声明: 本站内容部分来源网络,版权归作者所有,如有侵权,请联系我们删除!