python - 使用python的三对角粗麻布

Question

我正在研究 Python 中数值三对角粗麻布的实现。为此，我使用有限差分法，这是我迄今为止尝试过的：

import numpy as np

def principal_diagonal(x, fun):
    h = 10e-6
    result = list()
    for i in range(len(x)):
        
        temp_forward = np.array(x)
        temp_backward = np.array(x)
        temp_forward[i] = temp_forward[i] + h
        temp_backward[i] = temp_backward[i] - h
        value = (fun(temp_forward) - 2 * fun(x) + fun(temp_backward)) / (h ** 2)
        result.append(value)
    return(np.array(result))

def one_above_diagonal(x, fun):
    h =  10e-6 
    result = list()
    for i in range(len(x)-1):
        clean = np.array(x)
        forward_first = np.array(x)
        forward_first[i] = forward_first[i] + h
        forward_second = np.array(forward_first)
        forward_second[i+1] = forward_second[i+1] + h
        clean[i+1] = clean[i+1] + h
        value =  ((fun(forward_second) - fun(clean)) - (fun(forward_first) - fun(x))) / (h*h)
        result.append(value)
    return(np.array(result))

def tridiag_hessian(x, fun):
    new_array = np.ndarray((len(x), len(x)))
    di = np.diag_indices(len(x))
    new_array[di] = principal_diagonal(x, fun)
    new_array[(np.arange(len(x)-1), np.arange(1,len(x)))] = one_above_diagonal(x, fun)
    new_array[(np.arange(1,len(x)), np.arange(len(x)-1))] = one_above_diagonal(x, fun)
    return(new_array)

我已经使用 Rosenbrock 函数测试了我的函数，结果如下：

from scipy.optimize import rosen
tridiag_hessian(np.array([1, 2], "float64"), rosen)

array([[ 401.99964246, -400.00216472],
       [-400.00216472,  200.00001655]])

结果与实际解决方案非常接近，这是不错的，但对于更复杂的问题，我希望有更高的精度。我的代码中有什么可以调整的吗？或者，也许，另一种方法？我尝试过使用numdifftools然后继续对该函数的结果进行“三对角化”，但是当我在涉及具有超过 40 个变量的函数的问题中尝试这种方法时，速度非常慢。我也考虑过使用np.sqrt(np.finfo(float).eps)而不是 0.000001，但我没有获得任何精度。

标签： pythonnumpynumerical-methodsnumerical-analysis