Thomas求解算法Python 3

Question

我需要编写一个实现Thomas算法的程序。我的程序编译得很好，但是当我比较np.linalg.求解numpy im中的函数对同一个系统得到不同的结果

这是我的密码

import numpy as np

## Tri Diagonal Matrix Algorithm(a.k.a Thomas algorithm) solver
def TDMAsolver(a, b, c, d):

    nf = len(d) # number of equations
    ac, bc, cc, dc = map(np.array, (a, b, c, d)) # copy arrays
    for it in range(1, nf):
        mc = ac[it-1]/bc[it-1]
        bc[it] = bc[it] - mc*cc[it-1] 
        dc[it] = dc[it] - mc*dc[it-1]

    xc = bc
    xc[-1] = dc[-1]/bc[-1]

    for il in range(nf-2, -1, -1):
        xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il]

    return xc

下面是我如何测试它

^{pr2}$

输出是

[5 1 0 2]

当我这么做的时候

a = np.array([3,3,3,3])
c = np.array([2,2,2,2])
d = np.array([2,2,2,2])
b = np.array([12,17,14,7])

A = np.diag(a,0) + np.diag(c[1:],-1) + np.diag(d[1:],1)
x = np.linalg.solve(A,b)
print(x)

输出是

[ 2.  3.  2.  1.]

你猜为什么？在

谢谢你

Answer 1

np.linalg.solve调用的结果就是您需要得到的结果。在

您的功能有两个问题：

• 当使用int类型的numpy数组调用它时，对（它们的副本）的所有赋值也将导致int，这意味着在中间计算中会丢失小数部分。

• 它使用它的第二个参数作为主对角线（即[2,2,2,2]），而传递给np.linalg.solve的矩阵，其中[3,3,3,3]是主对角线。

我还建议在全局范围内使用与函数中相同的命名约定，因为现在有了它，在全局范围内调用方程的右边b，但在函数中它被称为d。全局范围内的d是函数中的b。这无助于获取代码的夹点。在

1。避免截断为整数

当你看这个作业时，第一个问题就变得很明显了：

bc[it] = bc[it] - mc*cc[it-1] 

在第一次迭代中，计算结果为

^{pr2}$

它应该是1.66666667，但是如果在赋值之后检查bc[it]的值，就会看到它的值是1。这是因为numpy数组包含整型值。在

这可以通过替换以下内容来解决：

ac, bc, cc, dc = map(np.array, (a, b, c, d)) 

有：

ac, bc, cc, dc = (x.astype(float) for x in (a, b, c, d))

2。命名并传递变量right

要求解的方程如下所示：

( 3 2 0 0 )   ( x0 )   ( 12 )
( 2 3 2 0 ) . ( x1 ) = ( 17 )
( 0 2 3 2 )   ( x2 )   ( 14 )
( 0 0 2 3 )   ( x3 )   (  7 )

您的函数期望这些值出现在参数中，如下所示

( b0 c0 0  0  )   ( x0 )   ( d0 )
( a0 b1 c1 0  ) . ( x1 ) = ( d1 )
( 0  a1 b2 c2 )   ( x2 )   ( d2 )
( 0  0  a2 b3 )   ( x3 )   ( d3 )

有两个数组的大小应该与数组的大小相协调：

a = np.array([2,2,2])
b = np.array([3,3,3,3]) # main diagonal
c = np.array([2,2,2])
d = np.array([12,17,14,7]) # right side of equation

x = TDMAsolver(a,b,c,d) # pass in same order

print(x)

为了与命名约定和更短的数组相匹配，请对代码的验证部分进行如下调整：

A = np.diag(b,0) + np.diag(a,-1) + np.diag(c,1)
x = np.linalg.solve(A,b)
print(x)

这两种情况下的解决方案是：

x == [ 2.  3.  2.  1.]

看到它运行在ideone.com

Answer 2

似乎您正纠结于对角线和右列数组目的地（命名、赋值？）。在

如果b是主对角线，a是下对角线，c是上对角线，d是右栏

b c 0 0       d
a b c 0       d 
0 a b c       d
0 0 a b       d

那么你的代码看起来是对的。但是矩阵A的收缩是为了什么呢利纳格。解决？你能打印出现成的矩阵并检查元素位置吗？在