Почему эта реализация алгоритма Block Tridiagonal Thomas дает такие большие ошибки? ⇐ Python

[Anonymous]

См. ниже мою попытку реализации блочного трехдиагонального алгоритма Томаса. Однако если вы запустите это, вы получите относительно большую (10^-2) ошибку в блоке TMDA по сравнению с прямым решением np (10^-15), даже для этого очень простого случая. Более сложные тестовые примеры дают большие ошибки - думаю, ошибки начинают расти при обратной подстановке. Любая помощь относительно того, почему будет очень признательна!
import numpy as np
import torch

def solve_block_tridiagonal(a, b, c, d):

N = len(b)
x = np.zeros_like(d)

# Forward elimination with explicit C* and d* storage
C_star = np.zeros_like(c)
d_star = np.zeros_like(d)

# Initial calculations for C_0* and d_0*
C_star[0] = np.linalg.solve(b[0], c[0])
d_star[0] = np.linalg.solve(b[0], d[0])

# Forward elimination
for i in range(1, N - 1):
C_star = np.linalg.solve(b - a[i-1] @ C_star[i-1], c)
d_star = np.linalg.solve(b - a[i-1] @ C_star[i-1], d - a[i-1] @ d_star[i-1])

# Last d_star update for the last block
d_star[-1] = np.linalg.solve(b[-1] - a[-2] @ C_star[-2], d[-1] - a[-2] @ d_star[-2])

# Backward substitution
x[-1] = d_star[-1]
for i in range(N-2, -1, -1):
x = d_star - C_star @ x[i+1]

return x

def test_block_tridiagonal_solver():

N = 4

a = np.array([
[[1, 0.5], [0.5, 1]],
[[1, 0.5], [0.5, 1]],
[[1, 0.5], [0.5, 1]]
], dtype=np.float64)

b = np.array([
[[5, 0.5], [0.5, 5]],
[[5, 0.5], [0.5, 5]],
[[5, 0.5], [0.5, 5]],
[[5, 0.5], [0.5, 5]]
], dtype=np.float64)

c = np.array([
[[1, 0.5], [0.5, 1]],
[[1, 0.5], [0.5, 1]],
[[1, 0.5], [0.5, 1]]
], dtype=np.float64)

d = np.array([
[1, 2],
[2, 3],
[3, 4],
[4, 5]
], dtype=np.float64)

x = solve_block_tridiagonal(a, b, c, d)

# Construct the equivalent full matrix A_full and right-hand side d_full
A_full = np.block([
[b[0], c[0], np.zeros((2, 2)), np.zeros((2, 2))],
[a[0], b[1], c[1], np.zeros((2, 2))],
[np.zeros((2, 2)), a[1], b[2], c[2]],
[np.zeros((2, 2)), np.zeros((2, 2)), a[2], b[3]]
])

d_full = d.flatten() # Flatten d for compatibility with the full system

# Solve using numpy's direct solve for comparison
x_np = np.linalg.solve(A_full, d_full).reshape((N, 2))
# Print the solutions for comparison
print("Solution x from block tridiagonal solver (TMDA):\n", x, "\nResidual:", torch.sum(torch.abs(torch.tensor(A_full)@torch.tensor(x).flatten() - torch.tensor(d).flatten())))
print("Solution x from direct full matrix solver:\n", x_np, "\nResidual np:", torch.sum(torch.abs(torch.tensor(A_full)@torch.tensor(x_np).flatten() - torch.tensor(d).flatten())))
# Run the test function
test_block_tridiagonal_solver()


Подробнее здесь: https://stackoverflow.com/questions/791 ... -give-such