圆月山庄资源网

python实现三次样条插值
脚本专栏 2024/11/6 佚名
3 2
圆月山庄资源网 Design By www.vgjia.com

本文实例为大家分享了python实现三次样条插值的具体代码，供大家参考，具体内容如下

函数：

python实现三次样条插值

算法分析

三次样条插值。就是在分段插值的一种情况。

要求：

  • 在每个分段区间上是三次多项式（这就是三次样条中的三次的来源）
  • 在整个区间（开区间）上二阶导数连续（当然啦，这里主要是强调在节点上的连续）
  • 加上边界条件。边界条件只需要给出两个方程。构建一个方程组，就可以解出所有的参数。

这里话，根据第一类样条作为边界。（就是知道两端节点的导数数值，然后来做三次样条插值）

但是这里也分为两种情况，分别是这个数值是随便给的一个数，还是说根据函数的在对应点上数值给出。

情况一：两边导数数值给出

这里假设数值均为1。即 f′(x0)=f′(xn)=f′(xn)=1的情况。

情况一图像

python实现三次样条插值

情况一代码

import numpy as np
from sympy import *
import matplotlib.pyplot as plt


def 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 I


def ff(x): # f[x0, x1, ..., xk]
 ans = 0
 for i in range(len(x)):
 temp = 1
 for j in range(len(x)):
  if i != j:
  temp *= (x[i] - x[j])
 ans += f(x[i]) / temp
 return ans


def calM():
 lam = [1] + [1 / 2] * 9
 miu = [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) * 2
 for 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.I
 return Ms.tolist()[0]


def calnf(x):
 nf = []
 for i in range(len(x) - 1):
 nf.append(cal(x[i], x[i + 1], i))
 return nf


def calf(f, x):
 y = []
 for i in x:
 y.append(f.subs(n, i))
 return y


def nfSub(x, nf):
 tempx = np.array(range(11)) - 5
 dx = []
 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'] = False
 x = np.linspace(-5, 5, 101)
 y = f(x)
 Ly = nfSub(x, nf)
 plt.plot(x, y, label='原函数')
 plt.plot(x, Ly, label='三次样条插值函数')
 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 - y
 temp = abs(temp)
 print(temp.mean())


if __name__ == '__main__':
 x = np.array(range(11)) - 5
 y = f(x)

 n, m = symbols('n m')
 init_printing(use_unicode=True)
 Ms = calM()
 nf = calnf(x)
 draw(nf)
 lossCal(nf)

情况二：两边导数数值由函数本身算出

这里假设数值均为1。即 f′(xi)=S′(xi)(i=0,n)f′(xi)=S′(xi)(i=0,n)的情况。

情况二图像

python实现三次样条插值

情况二代码

import numpy as np
from sympy import *
import matplotlib.pyplot as plt


def 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 I


def ff(x): # f[x0, x1, ..., xk]
 ans = 0
 for i in range(len(x)):
 temp = 1
 for j in range(len(x)):
  if i != j:
  temp *= (x[i] - x[j])
 ans += f(x[i]) / temp
 return ans


def calM():
 lam = [1] + [1 / 2] * 9
 miu = [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) * 2
 for 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.I
 return Ms.tolist()[0]


def calnf(x):
 nf = []
 for i in range(len(x) - 1):
 nf.append(cal(x[i], x[i + 1], i))
 return nf


def calf(f, x):
 y = []
 for i in x:
 y.append(f.subs(n, i))
 return y


def nfSub(x, nf):
 tempx = np.array(range(11)) - 5
 dx = []
 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'] = False
 x = np.linspace(-5, 5, 101)
 y = f(x)
 Ly = nfSub(x, nf)
 plt.plot(x, y, label='原函数')
 plt.plot(x, Ly, label='三次样条插值函数')
 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 - y
 temp = abs(temp)
 print(temp.mean())


if __name__ == '__main__':
 x = np.array(range(11)) - 5
 y = f(x)

 n, m = symbols('n m')
 init_printing(use_unicode=True)
 Ms = calM()
 nf = calnf(x)
 draw(nf)
 lossCal(nf)

以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持。

标签：
python,三次样条插值,分段插值
圆月山庄资源网 Design By www.vgjia.com
广告合作：本站广告合作请联系QQ：858582 申请时备注：广告合作（否则不回）
免责声明：本站文章均来自网站采集或用户投稿，网站不提供任何软件下载或自行开发的软件！ 如有用户或公司发现本站内容信息存在侵权行为，请邮件告知！ 858582#qq.com
圆月山庄资源网 Design By www.vgjia.com
www.vgjia.com 圆月山庄资源网
139,976互联网资源
144,792高清电影
21,817无损音乐
631,128技术资源

更新日志

友情链接

杰晶网络 DDR爱好者之家 桃源资源网 杰网资源 富贵资源网 南强小屋 铁雪资源网 幽灵资源网 万梅资源网 狼山资源网 白云岛资源网 昆仑资源网 相思资源网 明霞山资源网 内蒙古资源网 黑松山资源网 茶园资源网 饿虎岗资源网 大旗谷资源网 常春岛资源网 岱庙资源网 兴国资源网 快活林资源网 蝙蝠岛资源网 帝王谷资源网 白云城资源网 伏龙阁资源网 清风细雨楼 天枫庄资源网 圆月山庄资源网 无争山庄资源网 神水资源网 移花宫资源网 神剑山庄资源网 无为清净楼资源网 金钱帮资源网 丐帮资源网 华山资源网 极乐门资源网 小李飞刀资源网 凤求凰客栈 风云阁资源网 金狮镖局 鸳鸯亭资源网 千金楼资源网 更多链接