Cubic Spline Код Python, производящий линейные сплайны

Question

Редактировать: Я не ищу вас для отладки этого кода. Если вы знакомы с этим известным алгоритмом, то вы можете помочь. Обратите внимание, что алгоритм корректно генерирует коэффициенты.

Этот код для кубической сплайновой интерполяции производит линейные сплайны, и я не могу понять, почему (пока). Алгоритм исходит из численного анализа Бурдена, который примерно идентичен псевдокоду 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)) 

Для скучающих или сверхуспевающих ...

Вот код для тестирования, интервал х: [1, 9], у: [0, 19,7750212]. Функция тест XLN (х), поэтому мы начинаем 1 и увеличение на .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,6х^(- 2x) sin (3 * пи * х) между 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 

---

• Math, вероятно, следует сделать с помощью numpy или scipy.

---

Кроме того, you should read Idiomatic Python by David Goodger.