큐빅 스플라인 선형 스플라인을 생성하는 파이썬 코드

Question

편집 :이 코드를 디버깅하지 않으려 고합니다. 이 잘 알려진 알고리즘에 익숙하다면 도움이 될 것입니다. 알고리즘은 계수를 올바르게 생성합니다.

이 입방 스플라인 보간 코드는 선형 스플라인을 생성하므로 왜 (아직) 그럴 수없는 것입니까? 이 알고리즘은 Burden의 수치 해석 (Numerical Analysis)에서 나온 것으로 의사 코드 here과 거의 동일하거나 주석의 링크에서 해당 책을 찾을 수 있습니다 (3 장 참조). 코드가 올바른 계수를 생성하고 있습니다. 구현을 오해하고 있다고 생각합니다. 모든 의견을 크게 주시면 감사하겠습니다. 또한, 프로그래밍에 익숙하지 않아 코딩이 얼마나 나쁜지에 대한 피드백도 환영합니다. 나는 h, a, c의 관점에서 선형 시스템의 사진을 업로드하려했지만 새로운 사용자로서 나는 할 수 없다. 알고리즘이 해결하고 var alpha에 의해 설정되는 3 중선 선형 시스템의 시각을 원한다면 책의 주석에있는 링크를 참조하십시오 (3 장 참조). 시스템은 대각선으로 지배적이므로 시스템을 잘 알고 있습니다. 고유 한 해 c0, ..., cn이 존재한다. ci 값을 알면 다른 계수들이 따라옵니다. 뚫린 또는 overachieving 들어

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, 10], Y [0 19.7750212]. 테스트 기능은 XLN (x)를, 그래서 우리는 ... 그리고 플로팅 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 

까지 0.1로 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 %)의 계수를 먼저 구성한 다음이를 구현하는 방법으로 완전히 기능하는 3 차 스플라인 보간입니다. 분명히 이것이 유일한 방법은 아닙니다. 나는 다른 접근법에 대해 연구하고 관심이 있다면 게시 할 수 있습니다. 코드를 명확하게하는 한 가지 방법은 해결되는 선형 시스템의 이미지가 될 것이지만 내 담당자가 최대 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^(- 2 배) 0과 1 사이의 죄 (3 * 파이 *의 X) :

# 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): 
      """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.""" 
  

  ---

  • 당신은 반복 요소의 목록을 만들 수 있습니다
    이

    def cubic_spline(xx,yy): 
    

    같은를 작성 해주세요 : 대신

  . 이처럼

은 :

>>> [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.