0

我正在研究函数逼近，在尝试理解/实现多项式插值时，我在这里找到了一个示例。我发现下面的代码是一个很好的例子来理解实际发生的事情，而不是使用准备好的函数，但是它没有运行：

定义插值算法。本质上，我们试图将真实 f 表示为基函数 (psi-s) 的线性组合。

    import sympy as sym
    def interpolation(f, psi, points):
    N = len(psi) - 1 #order of the approximant polynomial
    A = sym.zeros((N+1, N+1)) # initiating the square matrix, whose regular element is psi evaluated at each nodes
    b = sym.zeros((N+1, 1)) # original function f evaluated at the selected nodes
    psi_sym = psi  # save symbolic expression
    # Turn psi and f into Python functions
    x = sym.Symbol('x')
    psi = [sym.lambdify([x], psi[i]) for i in range(N+1)]
    f = sym.lambdify([x], f)
    for i in range(N+1):
        for j in range(N+1):
            A[i,j] = psi[j](points[i])
        b[i,0] = f(points[i])
    c = A.LUsolve(b) #finding the accurate weights for each psi
    # c is a sympy Matrix object, turn to list
    c = [sym.simplify(c[i,0]) for i in range(c.shape[0])]
    u = sym.simplify(sum(c[i,0]*psi_sym[i] for i in range(N+1)))
    return u, c

真函数 f:= 10(x-1)^2 -1，节点：x0:= 1 + 1/3 和 x1 = 1 + 2/3。间隔：[1,2]。

x = sym.Symbol('x')
f = 10*(x-1)**2 - 1
psi = [1, x] # approximant polynomial of order 1 (linear approximation)
Omega = [1, 2] #interval
points = [1 + sym.Rational(1,3), 1 + sym.Rational(2,3)]
u, c = interpolation(f, psi, points)
comparison_plot(f, u, Omega)

代码不运行。错误发生在行

A = sym.zeros((N+1, N+1))

错误消息：ValueError: (2, 2) is not an integer
但是 A 不应该是整数，它是一个方阵，其每个元素在每个节点上都经过 psi 评估。f = A*c。

谢谢！！！

4

0 回答 0