函數(shù)

$$y=\frac{1}{1+x^2}$$

算法分析

三次樣條插值。就是在分段插值的一種情況。
要求：

• 在每個分段區(qū)間上是三次多項式（這就是三次樣條中的三次的來源）
• 在整個區(qū)間（開區(qū)間）上二階導(dǎo)數(shù)連續(xù)（當(dāng)然啦，這里主要是強(qiáng)調(diào)在節(jié)點上的連續(xù)）
• 加上邊界條件。邊界條件只需要給出兩個方程。構(gòu)建一個方程組，就可以解出所有的參數(shù)。

這里話，根據(jù)第一類樣條作為邊界。（就是知道兩端節(jié)點的導(dǎo)數(shù)數(shù)值，然后來做三次樣條插值）

但是這里也分為兩種情況，分別是這個數(shù)值是隨便給的一個數(shù)，還是說根據(jù)函數(shù)的在對應(yīng)點上數(shù)值給出。

情況一：兩邊導(dǎo)數(shù)數(shù)值給出

這里假設(shè)數(shù)值均為1。即 $f^{\prime }(x_0)=f^{\prime }(x_n)=1$的情況。

情況一圖像

情況一代碼

import numpy as np from sympy import * import matplotlib.pyplot as pltdef f(x):return 1 / (1 + x ** 2)def cal(begin, end, i):by = f(begin)ey = f(end)I = Ms[i] * ((end - n) ** 3) / 6 + Ms[i + 1] * ((n - begin) ** 3) / 6 + (by - Ms[i] / 6) * (end - n) + (ey - Ms[i + 1] / 6) * (n - begin)return Idef ff(x): # f[x0, x1, ..., xk]ans = 0for i in range(len(x)):temp = 1for j in range(len(x)):if i != j:temp *= (x[i] - x[j])ans += f(x[i]) / tempreturn ansdef calM():lam = [1] + [1 / 2] * 9miu = [1 / 2] * 9 + [1]# Y = 1 / (1 + n ** 2)# df = diff(Y, n)x = np.array(range(11)) - 5# ds = [6 * (ff(x[0:2]) - df.subs(n, x[0]))]ds = [6 * (ff(x[0:2]) - 1)]for i in range(9):ds.append(6 * ff(x[i: i + 3]))# ds.append(6 * (df.subs(n, x[10]) - ff(x[-2:])))ds.append(6 * (1 - ff(x[-2:])))Mat = np.eye(11, 11) * 2for i in range(11):if i == 0:Mat[i][1] = lam[i]elif i == 10:Mat[i][9] = miu[i - 1]else:Mat[i][i - 1] = miu[i - 1]Mat[i][i + 1] = lam[i]ds = np.mat(ds)Mat = np.mat(Mat)Ms = ds * Mat.Ireturn Ms.tolist()[0]def calnf(x):nf = []for i in range(len(x) - 1):nf.append(cal(x[i], x[i + 1], i))return nfdef calf(f, x):y = []for i in x:y.append(f.subs(n, i))return ydef nfSub(x, nf):tempx = np.array(range(11)) - 5dx = []for i in range(10):labelx = []for j in range(len(x)):if x[j] >= tempx[i] and x[j] < tempx[i + 1]:labelx.append(x[j])elif i == 9 and x[j] >= tempx[i] and x[j] <= tempx[i + 1]:labelx.append(x[j])dx = dx + calf(nf[i], labelx)return np.array(dx)def draw(nf):plt.rcParams['font.sans-serif'] = ['SimHei']plt.rcParams['axes.unicode_minus'] = Falsex = np.linspace(-5, 5, 101)y = f(x)Ly = nfSub(x, nf)plt.plot(x, y, label='原函數(shù)')plt.plot(x, Ly, label='三次樣條插值函數(shù)')plt.xlabel('x')plt.ylabel('y')plt.legend()plt.savefig('1.png')plt.show()def lossCal(nf):x = np.linspace(-5, 5, 101)y = f(x)Ly = nfSub(x, nf)Ly = np.array(Ly)temp = Ly - ytemp = abs(temp)print(temp.mean())if __name__ == '__main__':x = np.array(range(11)) - 5y = f(x)n, m = symbols('n m')init_printing(use_unicode=True)Ms = calM()nf = calnf(x)draw(nf)lossCal(nf)

情況二：兩邊導(dǎo)數(shù)數(shù)值由函數(shù)本身算出

這里假設(shè)數(shù)值均為1。即 $f^{\prime }(x_i)=S^{\prime }(x_i)(i=0,n)$的情況。

情況二圖像

情況二代碼

import numpy as np from sympy import * import matplotlib.pyplot as pltdef f(x):return 1 / (1 + x ** 2)def cal(begin, end, i):by = f(begin)ey = f(end)I = Ms[i] * ((end - n) ** 3) / 6 + Ms[i + 1] * ((n - begin) ** 3) / 6 + (by - Ms[i] / 6) * (end - n) + (ey - Ms[i + 1] / 6) * (n - begin)return Idef ff(x): # f[x0, x1, ..., xk]ans = 0for i in range(len(x)):temp = 1for j in range(len(x)):if i != j:temp *= (x[i] - x[j])ans += f(x[i]) / tempreturn ansdef calM():lam = [1] + [1 / 2] * 9miu = [1 / 2] * 9 + [1]Y = 1 / (1 + n ** 2)df = diff(Y, n)x = np.array(range(11)) - 5ds = [6 * (ff(x[0:2]) - df.subs(n, x[0]))]# ds = [6 * (ff(x[0:2]) - 1)]for i in range(9):ds.append(6 * ff(x[i: i + 3]))ds.append(6 * (df.subs(n, x[10]) - ff(x[-2:])))# ds.append(6 * (1 - ff(x[-2:])))Mat = np.eye(11, 11) * 2for i in range(11):if i == 0:Mat[i][1] = lam[i]elif i == 10:Mat[i][9] = miu[i - 1]else:Mat[i][i - 1] = miu[i - 1]Mat[i][i + 1] = lam[i]ds = np.mat(ds)Mat = np.mat(Mat)Ms = ds * Mat.Ireturn Ms.tolist()[0]def calnf(x):nf = []for i in range(len(x) - 1):nf.append(cal(x[i], x[i + 1], i))return nfdef calf(f, x):y = []for i in x:y.append(f.subs(n, i))return ydef nfSub(x, nf):tempx = np.array(range(11)) - 5dx = []for i in range(10):labelx = []for j in range(len(x)):if x[j] >= tempx[i] and x[j] < tempx[i + 1]:labelx.append(x[j])elif i == 9 and x[j] >= tempx[i] and x[j] <= tempx[i + 1]:labelx.append(x[j])dx = dx + calf(nf[i], labelx)return np.array(dx)def draw(nf):plt.rcParams['font.sans-serif'] = ['SimHei']plt.rcParams['axes.unicode_minus'] = Falsex = np.linspace(-5, 5, 101)y = f(x)Ly = nfSub(x, nf)plt.plot(x, y, label='原函數(shù)')plt.plot(x, Ly, label='三次樣條插值函數(shù)')plt.xlabel('x')plt.ylabel('y')plt.legend()plt.savefig('1.png')plt.show()def lossCal(nf):x = np.linspace(-5, 5, 101)y = f(x)Ly = nfSub(x, nf)Ly = np.array(Ly)temp = Ly - ytemp = abs(temp)print(temp.mean())if __name__ == '__main__':x = np.array(range(11)) - 5y = f(x)n, m = symbols('n m')init_printing(use_unicode=True)Ms = calM()nf = calnf(x)draw(nf)lossCal(nf)

總結(jié)

以上是生活随笔為你收集整理的三次样条插值Python实现的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。