edit: Я не ищу, чтобы вы отлаживали этот код. Если вы знакомы с этим хорошо известным алгоритмом, то вы можете помочь. Обратите внимание, что алгоритм правильно выдает коэффициенты.
Этот код для кубической интерполяции сплайнов создает линейные сплайны, и я не могу понять, почему (пока). Алгоритм основан на численном анализе Бурдена, который примерно идентичен псевдокоду здесь , или вы можете найти эту книгу по ссылке в комментариях (см. Главу 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()