病态线性方程组¶
约 612 个字 73 行代码 预计阅读时间 3 分钟
对于大型线性方程组，在计算中不断有误差累积，造成解丢失两、三位有效数字的情况是时有发生的。对于这种 丢失精度 的情况，一个简单可靠的方法可以使解能迅速恢复到计算机精度，即采用解的 迭代改进 方法，也可用于求解*病态线性方程组*。
计算原理¶
设线性方程组
$$ \mathbf{A x} = \mathbf{b} $$
的精确解为 \(\mathbf{x}\)，但实际计算中只能得到带误差的解 \(\mathbf{x} + \delta \mathbf{x}\)。
将误差解代入原方程会导致右向量 \(\mathbf{b}\) 产生摄动：
$$ \mathbf{A}(\mathbf{x} + \delta \mathbf{x}) = \mathbf{b} + \delta \mathbf{b} $$
由此可建立 残差方程 ：
$$ \mathbf{A} \delta \mathbf{x} = \mathbf{A}(\mathbf{x} + \delta \mathbf{x}) - \mathbf{b} $$
迭代改进过程：
- 初始化阶段：使用高斯消元法计算初始近似解
- 迭代改进循环（默认最多5次）
- 计算当前解的残差向量：\(\delta \mathbf{b} = Ax^{(k)} - b\)
- 求解修正量方程：\(AΔx = \delta \mathbf{b}\)
- 更新解向量：\(x^{(k+1)} = x^{(k)} - Δx\)
- 终止条件
- 达到预设的最大迭代次数
虽然无法求得完全精确的改正量 \(\delta \mathbf{x}\)，但经过一次迭代后改进量会显著减小。实践中通常循环改进1-2次即可获得理想结果。
Fortran数值实现¶
subroutine iterprove(a,b,x,n)
! 对线性方程组Ax=b的直接法进行迭代改进
! 计算过程：
! 1. 调用选主元消去法计算
! 2. 计算残差方程
! 3. 迭代改进
integer, intent(in) :: n
real*8, intent(in) :: a(n,n), b(n)
real*8, intent(out) :: x(n)
real*8 :: x1(n), x2(n), db(n), dx(n) ! db为残差向量
integer :: i
integer :: iter_max = 5 ! 默认迭代次
character(len=32) :: fmt_str
write(fmt_str, '(a,I0,a)') '(" ",I4,2X,', n, 'F10.6)'
!先调用一次计算，得到初值
call gauss(A,b,x1,N)
! x1=1d0 ! 若由gauss得到的初值非常差，可以手动设置初值
write(*,'(/,T25,"Iterative Improvement Sequence",//,"Iter",T10,"Solution Vector")')
write(*,fmt_str) 0, x1
do i=1,iter_max
db=matmul(A,x1)-b
call gauss(A,db,dx,N)
!由此得到改进量 dx
!x2 为改进后的新解
x2=x1-dx
!更新变量
x1=x2
write(*,fmt_str) i, x1
end do
x = x1
end subroutine iterprove
数值案例¶
主程序：
program main
use matrix
implicit none
integer, parameter :: n = 6
real*8 :: a(n,n), b(n), x(n)
! 使用reshape按行填充矩阵A和B
a = reshape([3.3435d0, 1.4946d0, 0.3834d0, -1.9537d0, -2.4013d0, -3.5446d0, &
1.0198d0, 6.1618d0, 4.9613d0, 2.7491d0, 3.0007d0, -3.6393d0, &
2.3703d0, 2.7102d0, 4.2182d0, 1.1730d0, -0.6859d0, 3.6929d0, &
1.5408d0, 4.1334d0, 0.5732d0, 5.6869d0, 3.1065d0, 0.7970d0, &
3.8921d0, 3.4762d0, 3.9335d0, -2.1556d0, -6.1815d0, 0.4986d0, &
2.4815d0, 8.2582d0, 4.6190d0, -0.0022d0, -0.3620d0, -8.5505d0], [n,n], order=[2,1])
b = [5.0537d0, 0.0228d0, -3.4177d0, 7.6380d0, 0.5575d0, 9.8252d0]
call iterprove(a, b, x, n)
end program main
输出：
Iterative Improvement Sequence
Iter Solution Vector
0 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
1 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
2 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
3 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
4 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
5 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
在这里是首次获取初值使用高斯求解，初值已经很好了，每次迭代，变化不是很明显。可以将初值设为1，再次迭代运行：
Iterative Improvement Sequence
Iter Solution Vector
0 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
2 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
3 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
4 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
5 0.461329 5.216169 -4.628165 -3.063036 1.374872 1.465094
稍加迭代，解向量就可以达到较为理想的效果。
这里使用的 迭代改进 与方程组的迭代求解技术是不同的概念。仍然采用直接方法计算线性方程，只是如果方程解不是太理想，可以通过求解残差方程以改进原先的计算结果，或者检验方程解的收敛性。
注意：子程序中默认：
implicit real*8(a-z)
参考文献¶
- 宋叶志,茅永兴,赵秀杰.Fortran 95/2003科学计算与工程[M].清华大学出版社,2011.