Python 三次样条插值法

代码

'''

本函数通过三次样条插值法进行函数值计算

'''

# 三次样条插值
import numpy as np

# 用于存放x,y,m的值
x = np.array([1,2,4,5])
y = np.array([1,3,4,2])
m = np.array([17/8,None,None,-19/8])
lens = len(x)
x_f = 3.0     # 待插值点
# 用于存放几种数据
h = np.zeros(lens-1)
a = np.zeros(lens)
b = np.zeros(lens)
n = np.zeros([lens,lens])

for i in range(len(h)):
    h[i] = x[i+1] - x[i]

for i in range(1,len(a)-1):
    a[i] = h[i-1]/(h[i]+h[i-1])
    b[i] = 3*(((1-a[i])*(y[i]-y[i-1]))/h[i-1]+(a[i]*(y[i+1]-y[i]))/h[i])
a[0] = 1
b[0] = (3*(y[1]-y[0]))/h[0]
b[-1] = (3*(y[-1]-y[-2]))/h[-1]

for i in range(len(n)):
    n[i][i] = 2
    if i != len(n) - 1:
        n[i][i + 1] = a[i]
    if i != 0:
        n[i][i-1] = 1 - a[i]

for i in range(len(n)-1):
    if i == 0:
        m[i+1] = (b[i] - n[i][i]*m[i])/a[i]
    else:
        m[i+1] = (b[i] - n[i][i]*m[i] - n[i][i-1]*m[i-1])/a[i]

p = 0
for i in range(len(x)-1):
    if x[i] < x_f <= x[i + 1]:
        p = i
        break

if p == 0:
    print('待求值点超出函数范围!')
    exit()

fun_1 = (1+2*(x_f-x[p])/(x[p+1]-x[p]))*(((x_f-x[p+1])/(x[p]-x[p+1]))**2)*y[p]
fun_2 = (1+2*(x_f-x[p+1])/(x[p]-x[p+1]))*(((x_f-x[p])/(x[p+1]-x[p]))**2)*y[p+1]
fun_3 = (x_f-x[p])*(((x_f-x[p+1])/(x[p]-x[p+1]))**2)*m[p]
fun_4 = (x_f-x[p+1])*(((x_f-x[p])/(x[p+1]-x[p]))**2)*m[p+1]

fun = fun_1 + fun_2 + fun_3 + fun_4

print('函数在点 {0} 处的插值函数值为: {1} '.format(x_f,fun))

评论 3 您还未登录,请先 登录 后发表或查看评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:1024 设计师:我叫白小胖 返回首页

打赏作者

努力写代码的瀟

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值