使用Taichi Linear solver求解方程结果与Numpy || matlab求解的不一致

太极开发者小伙伴,你们好,我目前遇到一个问题,我在求解一个线性方程Ku=F的时候,发现同样的K,F,使用taichi求解器计算出的u一直为0,与使用numpy或者matlab计算出的结果不一致。下面是我的代码和结果:

import taichi as ti
import numpy as np

Knumpy = np.array([    
        [     1,      0,      0,      0,      0,      0,      0,      0,      0,      0,      0,      0],
        [     0,      1,      0,      0,      0,      0,      0,      0,      0,      0,      0,      0],
        [     0,      0,  519.2,  115.4, -158.7, -108.2,      0,      0,  57.69,  57.69, -129.8, -79.33],
        [     0,      0,  115.4,  519.2, -7.212, -259.6,      0,      0,  57.69,  259.6, -79.33, -129.8],
        [     0,      0, -158.7, -7.212,  259.6,  7.212,      0,      0, -129.8,  21.63,  28.85, -21.63],
        [     0,      0, -108.2, -259.6,  7.212,  259.6,      0,      0,  21.63, -129.8,  79.33,  129.8],
        [     0,      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,      0],
        [     0,      0,      0,      0,      0,      0,      0,      1,      0,      0,      0,      0],
        [     0,      0,  57.69,  57.69, -129.8,  21.63,      0,      0,  519.2,  115.4, -158.7, -7.212],
        [     0,      0,  57.69,  259.6,  21.63, -129.8,      0,      0,  115.4,  519.2, -108.2, -259.6],
        [     0,      0, -129.8, -79.33,  28.85,  79.33,      0,      0, -158.7, -108.2,  259.6,  108.2],
        [     0,      0, -79.33, -129.8, -21.63,  129.8,      0,      0, -7.212, -259.6,  108.2,  259.6],
    ])

fnumpy = np.array([
    0,
    0,
    0,
    0,
    9.375,
    0,
    0,
    0,
    0,
    0,
    9.375,
    0
])

np.set_printoptions(precision=3, suppress=True)
np.set_printoptions(formatter={'float': '{: 0.3f}'.format})
print(np.linalg.inv(Knumpy) @ fnumpy)


ti.init(arch=ti.cpu)

f = ti.ndarray(ti.f64, 12)    
f[4]=9.375      
f[10]=9.375     


KBuilder = ti.linalg.SparseMatrixBuilder(12, 12, max_num_triplets=10000)

@ti.kernel
def buildSparseMatrix(K1: ti.types.sparse_matrix_builder(),K2: ti.types.ndarray()):
    for i, j in ti.ndrange(K2.shape[0],K2.shape[1]):
        K1[i,j] += K2[i,j]
@ti.kernel
def print_arr(arr : ti.types.ndarray(dtype=ti.f64, ndim=1)):
    for i in ti.grouped(arr):
        print(arr[i])


buildSparseMatrix(KBuilder,Knumpy)
kb = KBuilder.build()
solver = ti.linalg.SparseSolver(solver_type="LDLT")
solver.analyze_pattern(kb)
solver.factorize(kb)
dv = solver.solve(f)
print_arr(dv)

@houkensjtu 可以帮忙看下吗?

你好:wave: 这边的问题是需要显式指定一下 SparseMatrixSparseSolver 的数据精度:

KBuilder = ti.linalg.SparseMatrixBuilder(12, 12, max_num_triplets=10000, dtype=ti.f64)
solver = ti.linalg.SparseSolver(solver_type="LDLT", dtype=ti.f64)

以保证和你的方程右边精度一致。我确认了一下指定上述数据精度后结果应该是正确的。
出现这个问题的原因是 SparseMatrixSparseSolver 默认精度都是 f32,和你的右边精度不一致导致了求解出错了。