旋转行求解联立方程(针对特殊情况).(python,数值分析)

第六课 旋转行求解联立方程

在使用传统的高斯消元法或[L][U]因子分解法(程序2.2)求解非对称方程时，我们回避了系数矩阵的一个主对角分量应该从零开始还是在求解过程中变为零的问题。在下一个程序中，我们将说明如何通过使用通常称为行交换技术来解决这个问题。所以技术包括在尚未处理的行中搜索“最大”(绝对)系数，并通过行和列交换将其移动到对角位置。
例如，在这个等式中，第2行中的“最大”系数是-20，因此这与第1行互换，给出


进行第一步消元之后，变成

上面的最大系数为第三行第三列，所以交换为第二行第二列，为

得到最后的方程为

向后迭代法得到

这与之前的解决方案相同，但是未知量是以不同的顺序由前面的列互换决定的。这避免了主对角线中的第一个数为0的情况。
程序如下
一个主程序和两个子程序，子程序分别为求解顺序的子程序lupfac，和之后的’向后迭代法’的子程序lupsol
主程序

#使用行变换的LU分解
import numpy as np
import math
import B
n=3
row=np.zeros((n,1),dtype=np.int64)
sol=np.zeros((n,1))
a=np.array([[10,1,-5],[-20,3,20],[5,3,5]],dtype=np.float)
b=np.array([[1],[2],[6]],dtype=np.float)
print('系数矩阵')
print(a)
print('右手边向量',b)
B.lupfac(a,row)
B.lupsol(a,b,sol,row)
print('向后迭代顺序',row)
print('解向量',sol)   

lupfac

def lupfac(a,row):
  n=a.shape[0]
  for i in range(1,n+1):
    row[i-1,0]=i
  for i in range(1,n):
    ip=i
    pval=a[row[ip-1]-1,ip-1]
    for j in range(i+1,n+1):
      if abs(a[row[j-1]-1,i-1])>abs(pval):
        ip=j
        pval=a[row[j-1]-1,i-1]
    if abs(pval)<1.0e-20:
      print('方程解有奇异')
      exit
    ih=row[ip-1,0]
    row[ip-1,0]=row[i-1,0]
    row[i-1,0]=ih
    for j in range(i+1,n+1):
      ie= row[j-1,0]
      pivot=a[ie-1,i-1]/pval
      a[ie-1,i-1]=pivot
      irow=row[i-1,0]
      for k in range(i+1,n+1):
        a[ie-1,k-1]=a[ie-1,k-1]-a[irow-1,k-1]*pivot
  if abs(a[row[n-1]-1,n-1])<1.0e-10:
    print('方程有奇异解')
    exit

lupsol

def lupsol(a,b,sol,row):
  temp=np.zeros((a.shape[0],1))
  n=a.shape[0]
  temp[:,0]=b[:,0]
  for i in range(1,n+1):
    irow=row[i-1,0]
    total=b[irow-1,0]
    if i>1:
      for j in range(1,i):
        total=total-a[irow-1,j-1]*temp[row[j-1,0]-1,0]
      temp[irow-1,0]=total
  for i in range(n,0,-1):
    irow=row[i-1,0] 
    total=temp[irow-1,0]
    if i<n:
      for j in range(i+1,n+1):
        total=total-a[irow-1,j-1]*temp[row[j-1]-1,0]
    temp[irow-1,0]=total/a[irow-1,i-1]
  sol[:,0]=temp[row[:,0]-1,0]

终端输出结果如下：

程序结果与计算结果一致。