生成线性样条的三次样条Python代码

Question

编辑：我不希望您调试此代码。如果您熟悉这个著名的算法，那么您可能会有所帮助。请注意，算法生成的系数是正确的。

这个三次样条插值的代码正在生成线性样条，我似乎还不明白为什么。该算法来自于baund的数值分析，它与伪代码here几乎相同，或者您可以从注释中的链接找到那本书（参见第3章，无论如何值得拥有）。代码正在生成正确的系数；我认为我误解了实现。任何反馈都非常感谢。另外，我是编程新手，所以任何关于我的代码有多糟糕的反馈也是受欢迎的。我试着用h，a，c来上传线性系统的图片，但是作为一个新用户我不能。如果你想要一个算法求解的三对角线性系统的视觉效果，并且它是由var alpha设置的，请参阅本书注释中的链接，见第3章。这个系统是严格对角占优的，所以我们知道存在唯一的解c0，…，cn。一旦我们知道ci值，其他的系数就跟着来了。

import matplotlib.pyplot as plt

# need some zero vectors...
def zeroV(m):
    z = [0]*m
    return(z)

 #INPUT: n; x0, x1, ... ,xn; a0 = f(x0), a1 =f(x1), ... , an = f(xn).
 def cubic_spline(n, xn, a, xd):
    """function cubic_spline(n,xn, a, xd) interpolates between the knots
       specified by lists xn and a. The function computes the coefficients
       and outputs the ranges of the piecewise cubic splines."""        

    h = zeroV(n-1)

    # alpha will be values in a system of eq's that will allow us to solve for c
    # and then from there we can find b, d through substitution.
    alpha = zeroV(n-1)

    # l, u, z are used in the method for solving the linear system
    l = zeroV(n+1)
    u = zeroV(n)
    z = zeroV(n+1)

    # b, c, d will be the coefficients along with a.
    b = zeroV(n)     
    c = zeroV(n+1)
    d = zeroV(n)    

for i in range(n-1):
    # h[i] is used to satisfy the condition that 
    # Si+1(xi+l) = Si(xi+l) for each i = 0,..,n-1
    # i.e., the values at the knots are "doubled up"
    h[i] = xn[i+1]-xn[i]  

for i in range(1, n-1):
    # Sets up the linear system and allows us to find c.  Once we have 
    # c then b and d follow in terms of it.
    alpha[i] = (3./h[i])*(a[i+1]-a[i])-(3./h[i-1])*(a[i] - a[i-1])

# I, II, (part of) III Sets up and solves tridiagonal linear system...
# I   
l[0] = 1      
u[0] = 0      
z[0] = 0

# II
for i in range(1, n-1):
    l[i] = 2*(xn[i+1] - xn[i-1]) - h[i-1]*u[i-1]
    u[i] = h[i]/l[i]
    z[i] = (alpha[i] - h[i-1]*z[i-1])/l[i]

l[n] = 1
z[n] = 0
c[n] = 0

# III... also find b, d in terms of c.
for j in range(n-2, -1, -1):      
    c[j] = z[j] - u[j]*c[j+1]
    b[j] = (a[j+1] - a[j])/h[j] - h[j]*(c[j+1] + 2*c[j])/3.
    d[j] = (c[j+1] - c[j])/(3*h[j])   

# This is my only addition, which is returning values for Sj(x). The issue I'm having
# is related to this implemention, i suspect.
for j in range(n-1): 
    #OUTPUT:S(x)=Sj(x)= aj + bj(x - xj) + cj(x - xj)^2 + dj(x - xj)^3; xj <= x <= xj+1)
    return(a[j] + b[j]*(xd - xn[j]) + c[j]*((xd - xn[j])**2) + d[j]*((xd - xn[j])**3))

对于那些无聊的，或是有成就感的。。。

这是测试代码，间隔是x:[1，9]，y:[0，19.7750212]。测试函数是xln（x），所以我们从1开始，增加0.1到9。

ln = [] 
ln_dom = [] 
cub = [] 
step = 1. 
X=[1., 9.] 
FX=[0, 19.7750212]
while step <= 9.:
    ln.append(step*log(step))
    ln_dom.append(step)
    cub.append(cubic_spline(2, x, fx, step))
    step += 0.1

…对于密谋：

plt.plot(ln_dom, cub, color='blue')
plt.plot(ln_dom, ln, color='red')
plt.axis([1., 9., 0, 20], 'equal')
plt.axhline(y=0, color='black')
plt.axvline(x=0, color='black')
plt.show()

Answer 1

好的，开始工作了。问题出在我的实现上。我用了另一种方法，分别构造样条曲线，而不是连续构造。这是一种功能完备的三次样条插值方法，首先构造样条多项式的系数（这是99%的工作），然后实现它们。显然这不是唯一的办法。如果有兴趣的话，我可以用另一种方式发布。有一件事可以澄清代码是一个线性系统的图像，这是解决了，但我不能张贴图片，直到我的代表达到10。如果您想深入了解算法，请参阅上面注释中的“教科书”链接。

import matplotlib.pyplot as plt
from pylab import arange
from math import e
from math import pi
from math import sin
from math import cos
from numpy import poly1d

# need some zero vectors...
def zeroV(m):
    z = [0]*m
    return(z)

#INPUT: n; x0, x1, ... ,xn; a0 = f(x0), a1 =f(x1), ... , an = f(xn).
def cubic_spline(n, xn, a):
"""function cubic_spline(n,xn, a, xd) interpolates between the knots
   specified by lists xn and a. The function computes the coefficients
   and outputs the ranges of the piecewise cubic splines."""        

    h = zeroV(n-1)

    # alpha will be values in a system of eq's that will allow us to solve for c
    # and then from there we can find b, d through substitution.
    alpha = zeroV(n-1)

    # l, u, z are used in the method for solving the linear system
    l = zeroV(n+1)
    u = zeroV(n)
    z = zeroV(n+1)

    # b, c, d will be the coefficients along with a.
    b = zeroV(n)     
    c = zeroV(n+1)
    d = zeroV(n)    

    for i in range(n-1):
        # h[i] is used to satisfy the condition that 
        # Si+1(xi+l) = Si(xi+l) for each i = 0,..,n-1
        # i.e., the values at the knots are "doubled up"
        h[i] = xn[i+1]-xn[i]  

    for i in range(1, n-1):
        # Sets up the linear system and allows us to find c.  Once we have 
        # c then b and d follow in terms of it.
        alpha[i] = (3./h[i])*(a[i+1]-a[i])-(3./h[i-1])*(a[i] - a[i-1])

    # I, II, (part of) III Sets up and solves tridiagonal linear system...
    # I   
    l[0] = 1      
    u[0] = 0      
    z[0] = 0

    # II
    for i in range(1, n-1):
        l[i] = 2*(xn[i+1] - xn[i-1]) - h[i-1]*u[i-1]
        u[i] = h[i]/l[i]
        z[i] = (alpha[i] - h[i-1]*z[i-1])/l[i]

    l[n] = 1
    z[n] = 0
    c[n] = 0

    # III... also find b, d in terms of c.
    for j in range(n-2, -1, -1):      
        c[j] = z[j] - u[j]*c[j+1]
        b[j] = (a[j+1] - a[j])/h[j] - h[j]*(c[j+1] + 2*c[j])/3.
        d[j] = (c[j+1] - c[j])/(3*h[j]) 

    # Now that we have the coefficients it's just a matter of constructing
    # the appropriate polynomials and graphing.
    for j in range(n-1):
        cub_graph(a[j],b[j],c[j],d[j],xn[j],xn[j+1])

    plt.show()

def cub_graph(a,b,c,d, x_i, x_i_1):
    """cub_graph takes the i'th coefficient set along with the x[i] and x[i+1]'th
       data pts, and constructs the polynomial spline between the two data pts using
       the poly1d python object (which simply returns a polynomial with a given root."""

    # notice here that we are just building the cubic polynomial piece by piece
    root = poly1d(x_i,True)
    poly = 0
    poly = d*(root)**3
    poly = poly + c*(root)**2
    poly = poly + b*root
    poly = poly + a

    # Set up our domain between data points, and plot the function
    pts = arange(x_i,x_i_1, 0.001)
    plt.plot(pts, poly(pts), '-')
    return

如果您想测试它，这里有一些数据可以用来开始，这些数据来自 函数1.6e^（-2x）sin（3*pi*x）介于0和1之间：

# These are our data points
x_vals = [0, 1./6, 1./3, 1./2, 7./12, 2./3, 3./4, 5./6, 11./12, 1]

# Set up the domain
x_domain = arange(0,2, 1e-2)

fx = zeroV(10)

# Defines the function so we can get our fx values
def sine_func(x):
    return(1.6*e**(-2*x)*sin(3*pi*x))

for i in range(len(x_vals)):
    fx[i] = sine_func(x_vals[i])

# Run cubic_spline interpolant.
cubic_spline(10,x_vals,fx)

Answer 2

对编码风格的评论：

---

• 你的意见和文件在哪里？至少，提供函数文档，这样人们就可以知道您的函数应该如何使用。

而不是：

def cubic_spline(xx,yy):

请写下如下内容：

def cubic_spline(xx, yy):
    """function cubic_spline(xx,yy) interpolates between the knots
    specified by lists xx and yy. The function returns the coefficients
    and ranges of the piecewise cubic splines."""

---

• 通过在列表上使用*运算符，可以生成重复元素的列表。

像这样：

>>> [0] * 10
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

这样zeroV函数就可以被[0] * m替换。

不要用易变的类型！（尤其是列表）。

>>> inner_list = [1,2,3]
>>> outer_list = [inner_list] * 3
>>> outer_list
[[1, 2, 3], [1, 2, 3], [1, 2, 3]]
>>> inner_list[0] = 999
>>> outer_list
[[999, 2, 3], [999, 2, 3], [999, 2, 3]] # wut

---

• 数学计算应该使用numpy或scipy完成。

---

除此之外，you should read Idiomatic Python by David Goodger。