python三对角矩阵_追赶法求解三对角矩阵线性方程组的Python实现

追赶法是求大型稀疏方程之三对角线性方程组的三角分解方法，即求解方程组Ax=b，其中A为三对角矩阵，主对角线上的元素记为(a1, a2, …, an) ，紧贴主对角线上方的一根对角线上的元素记为(c1,…c(n-1) )，紧贴主对角线下方的一根对角线上的元素记为(d2, d3, …, dn)，b=(b1,…bn)

对三对角矩阵A进行Crout分解，有

A=[[α1 0 0 0 ....... [[1 β1 0 0 ......

γ1 α2 0 0 ...... * 0 1 β2 0 ...... = LU

0 γ1 α3 0 ......]] 0 0 1 β3......]]

依次求取αi，βi，其中

α1=a1，β1=c1/α1，γi=di

αi=ai-di*βi

βi=ci/αi

下面是一个例子

A=[[ 3 2 0

-1 3 2

0 -1 3...

..........]]的n阶三对角矩阵

b=(7,11,15,9)'

求解x

首先进行crout分解，将AX=b转化为LUx=b

令Ux=y，仿照平方根解法，利用Ly=b和Ux=y，解出x

代码如下

# -*- coding: utf-8 -*-

"""

Created on Mon Feb 11 20:19:00 2019

@author: 鹰皇

"""

#追赶法

import numpy as np

A=[[3.00000,2.00000,0.00000,0.00000],[-1.00000,3.00000,2.00000,0.00000],[0.00000,-1.00000,3.00000,2.00000],[0.00000,0.00000,-1.00000,3.00000]]

b=[[7.00000,11.00000,15.00000,9.00000]]

def get_base(A):#获得一个基，在上面修改得到答案

base=list(np.zeros((len(A),len(A))))

D=[]

for i in base:

D.append(list(i))

return D

def get_gamma(A)

4000

:#根据第一行公式γi=di

base=get_base(A)

for i in range(1,len(A)):

base[i][i-1]=A[i][i-1]

return base

def get_raw1(A,base):#根据第一行公式

base[0][0]=A[0][0]

base[0][1]=A[0][1]/base[0][0]

return base

def get_other_raws(A,base,i):#递推后面几行

base[i][i]=A[i][i]-A[i][i-1]*base[i-1][i]

base[i][i+1]=A[i][i+1]/base[i][i]

return base

def get_final_raw(A,base):#最后一行少一个β，另外求解，也可以和上面放在一起

base[-1][-1]=A[-1][-1]-A[-1][-2]*base[-2][-1]

return base

def get_all(A):#得到一个L和U并在一起的矩阵，由于U的主对角线为1，因此可以放在一起

base=get_base(A)

base=get_gamma(A)

base=get_raw1(A,base)

for i in range(1,len(A)-1):

get_other_raws(A,base,i)

get_final_raw(A,base)

return base

def get_lower(A):#获得L

for i in A:

for j in i:

if i.index(j)>A.index(i):

A[A.index(i)][i.index(j)]=0

return A

def get_upper(A):#获得U

for i in A:

for j in i:

if i.index(j)

A[A.index(i)][i.index(j)]=0.0

elif i.index(j)==A.index(i):

A[A.index(i)][i.index(j)]=1.0

return A

base=get_all(A)

lower=get_lower(base)

print np.mat(lower).I*np.mat(b).T#解得y

base=get_all(A)

upper=get_upper(base)

print np.mat(upper).I*np.mat(lower).I*np.mat(b).T解得x

最后解得x=(1,2,3,4)’，解答完成