Python实现第一型三次样条插值函数

例题

 

代码

import numpy as np
np.set_printoptions(suppress=True)

def h_u_lamb_A(x):
    h = []
    u = []
    lamb = []
    n=x.shape[0]
    A = np.mat(2*np.eye(n-2, n-2))
    A[0, 1] = 1
    A[n-3, n-4] = 1
    for i in range(0, n - 3):
        h.append(x[i + 2] - x[i+1])
    for i in range(0, n - 4):
        u.append(h[i] / (h[i] + h[i + 1]))
        lamb.append(h[i+1] / (h[i] + h[i + 1]))
        A[i + 1, i] = u[i]
        A[i + 1, i+2] = lamb[i]
    return h,u,lamb,A

def csb(x,f,j,dx0,dxn):
    n=x.shape[0]
    f0=np.zeros((j,x.shape[0]))
    if type(f) is np.ndarray:
        f0[0]=f.copy()
    else:
        for i in range(x.shape[0]):
            f0[0,i]=f(x[i])
    for i in range(1,j):
        for k in range(0,n-i):
            if i == 1 and k == i-1:
                f0[i, k] = dx0
            elif i == 1 and k == n-2:
                f0[i, k] = dxn
            else:
                f0[i,k]=(f0[i-1,k]-f0[i-1,k+1])/(x[k]-x[k+i])
    #print('所求%d阶差商表如下所示' % j,'\n',f0.T)
    #print('所求%d阶差商表如下所示' % j, '\n', f0.T[2][2])
    return f0.T

# e为误差
def Jacobi(A, b, x, e, times=100):
    length, width = np.shape(A)
    D = np.mat(np.diag(np.diag(A)))
    L = np.triu(A, 1)
    U = np.tril(A, -1)
    J = -D.I * (L + U)
    H = np.eye(length) - D.I * A
    eig, _ = np.linalg.eig(H)
    spectral_radius = max(abs(eig))
    if spectral_radius < 1:
        x0 = x
        x = J * x + D.I * b
        # x = D.I*(b-(L+U)*x0)
        k = 1
        while k < times:
            if abs(np.max(abs(x - x0), axis=0)) > e:
                x0 = x
                x = J * x + D.I * b
                k += 1
            else:   break
        return x

def Jacobi_resolve(A,x,csb_T):
    b = np.mat(np.zeros(x.shape[0] - 2)).I
    x0 = np.mat(np.zeros(x.shape[0] - 2)).I
    for i in range(0, x.shape[0] - 2):
        b[i, 0] = 6 * csb_T[i, 2]
    e = 0.001
    times = 100
    x0 = Jacobi(A, b, x0, e, times)
    return x0

def Calculate_Sx_k(xk,f,dx0,dxn):
    h, u, lamb, A = h_u_lamb_A(xk)
    csb_T = csb(xk, f, 3, dx0, dxn)  # 计算出查商表
    M = Jacobi_resolve(A, xk, csb_T)
    k = np.zeros((xk.shape[0] - 3, 4))  # 系数矩阵
    for j in range(0, xk.shape[0] - 3):
        k[j, 0] = f[j + 1]
        k[j, 1] = csb_T[j + 1, 1] - (M[j, 0] / 3 + M[j + 1, 0] / 6) * h[j]
        k[j, 2] = M[j, 0] / 2
        k[j, 3] = (M[j + 1, 0] - M[j, 0]) / (6 * h[j])
    return k

def Sx_desplay(xk,k):
    for j in range(0,xk.shape[0] - 3):
        print("S(x)=",'%.4f'%k[j, 0],"+",
              '%.4f'%k[j, 1],"* (x-" ,'%.4f'%xk[j + 1],")+",
              '%.4f'%k[j, 2],"* ((x-",'%.4f'%xk[j + 1],") ** 2) + ",
              '%.4f'%k[j, 3],"* ((x-",'%.4f'%xk[j + 1],") ** 3))  ",
              '%.4f'%xk[j + 1],"<=x<=",'%.4f'%xk[j + 2])

def Calculate_Sx(xk, k, x):
    for j in range(0, xk.shape[0] - 3):
        if x >= xk[j + 1] and x <= xk[j + 2]:
            Sx = '%.4f'%(k[j, 0] + k[j, 1] * (x - xk[j + 1]) + k[j, 2] * ((
                 x - xk[j + 1]) ** 2) + k[j, 3] * ((x - xk[j + 1]) ** 3))
    print("x=",x,"时,Sx=",Sx)


xk = np.array([0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 10])
f = np.array([2.51,2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80,5.80])
dx0 = 0.8
dxn=0.2
k = Calculate_Sx_k(xk,f,dx0,dxn)
Sx_desplay(xk,k)
for i in range(0,10):
    Calculate_Sx(xk, k, i + 0.5)