python插值程序_计算方法(2)——插值法(附Python程序)
给定一些数据,生成函数的方式有两种:插值,回归。
插值而得到的函数通过数据点,回归得到的函数不一定通过数据点。
下面给出拉格朗日插值,牛顿插值和Hermite插值的程序,
具体原理可参考课本,不再赘述。
拉格朗日插值法
线性插值 一次精度 需要2个节点
二次插值 二次精度 需要3个节点
n次插值 n次精度 需要n+1个节点
拉格朗日插值代码段(根据传入数据自动判断次数):
# 返回多项式
def p(x,a):
"""p(x,a)是x的函数,a是各幂次的系数"""
s = 0
for i in range(len(a)):
s += a[i]*x**i
return s
# n次拉格朗日插值
def lagrange_interpolate(x_list,y_list,x):
"""x_list 待插值的x元素列表y_list 待插值的y元素列表插值以后整个lagrange_interpolate是x的函数"""
if len(x_list) != len(y_list):
raise ValueError("list x and list y is not of equal length!")
# 系数矩阵
A = []
for i in range(len(x_list)):
A.append([])
for j in range(len(x_list)):
A[i].append(pow(x_list[i],j))
b = []
for i in range(len(x_list)):
b.append([y_list[i]])
# 求得各阶次的系数
a = lu_solve(A, b) # 用LU分解法解线性方程组,可以使用numpy的类似函数
a = transpose(a)[0] # change col vec a into 1 dimension
val = p(x,a)
print(x,val)
return val
其中lu_solve(A,b)是自己写的轮子,可以用numpy的numpy.linalg.sovle(A,b)来代替
到后面会有一期讲矩阵方程直接法,会有讲到如何写lu_solve()
看一看插值的效果如何
import numpy as np
import matplotlib.pyplot as plt
# 待插值的元素值
x_points = [0,1,2,3,4,5]
y_points = [1,5,4,8,7,12]
# 拉格朗日插值
x = np.linspace(0,5)
y = list(map(lambda t: lagrange_interpolate(x_points,y_points,t),x))
# 画图
plt.scatter(x_points,y_points,color = "orange")
plt.plot(x,y)
plt.legend(["lagrange interpolation","scattered points"])
plt.show()拉格朗日插值
可以看到,插值之后的函数可以较好地反映原数据点的情况。
牛顿插值法
涉及到的两个表,差分表和差商表:差分表
差商表
递归打印差分表
def difference_list(dlist): # Newton
if len(dlist)>0:
print(dlist)
prev,curr = 0,0
n = []
for i in dlist:
curr = i
n.append(curr - prev)
prev = i
n.pop(0)
difference_list(n)
打印差商表,并以列表的形式返回图中蓝色圈起来的部分
def difference_quotient_list(y_list,x_list = []):
if x_list == []:
x_list = [i for i in range(len(y_list))]
print(y_list)
prev_list = y_list
dq_list = []
dq_list.append(prev_list[0])
for t in range(1,len(y_list)):
prev,curr = 0,0
m = []
k = -1
for i in prev_list:
curr = i
m.append((curr - prev)/(x_list[k+t]-x_list[k]))
prev = i
k+=1
m.pop(0)
prev_list = m
dq_list.append(prev_list[0])
print(m)
return dq_list
牛顿插值代码段(调用了之前求差商表的代码)
def newton_interpolate(x_list,y_list,x):
coef = difference_quotient_list(y_list,x_list)
p = coef[0]
for i in range(1,len(coef)):
product = 1
for j in range(i):
product *= (x - x_list[j] )
p += coef[i]*product
return p
if __name__ == '__main__':
import numpy as np
import matplotlib.pyplot as plt
# 待插值的元素值
x_points = [0,1,2,3,4,5]
y_points = [1,5,4,8,7,12]
# 牛顿插值
x = np.linspace(0,5)
y = list(map(lambda t: newton_interpolate(x_points,y_points,t),x))
# 画图
plt.scatter(x_points,y_points,color = "orange")
plt.plot(x,y)
plt.legend(["newton interpolation","scattered points"])
plt.show()
插值效果图牛顿插值
可以看到,相同的插值节点,使用拉格朗日插值和牛顿插值的结果是相同的。
Hermite 插值法
三次Hermite插值不但要求在插值节点上插值多项式与函数值相等,还要求插值多项式的导数与函数导数相等。
进行Hermite插值需要六个信息:
这个比较好理解,通过
两点的插值有无限种方式。比如:
但是,通过限制住两个端点的导数值,就可以确定下来大致的插值形状。(对于Hermite插值,六个条件能确定出唯一的插值结果)
例如,可以编程求出
def hermite(x0,x1,y0,y1,y0_prime,y1_prime,x):
alpha0 = lambda x: ((x-x1)/(x0-x1))**2 * (2*(x-x0)/(x1-x0)+1)
alpha1 = lambda x: ((x-x0)/(x1-x0))**2 * (2*(x-x1)/(x0-x1)+1)
beta0 = lambda x: ((x-x1)/(x0-x1))**2 * (x-x0)
beta1 = lambda x: ((x-x0)/(x1-x0))**2 * (x-x1)
H = alpha0(x)*y0 + alpha1(x)*y1 + beta0(x)*y0_prime + beta1(x)*y1_prime
return H
if __name__ == '__main__':
import numpy as np
import matplotlib.pyplot as plt
f = lambda x: hermite(0,1,0,1,-1,-4,x)
x = np.linspace(0,1)
y = list(map(f,x))
plt.scatter([0,1],[0,1],color = "orange")
plt.plot(x,y)
plt.show()
龙格现象(Runge phenomenon)
然而,并不是所有的函数都可以通过提高插值次数来提高插值准确度。
比如著名的龙格函数
,从[-1,1]取10个点,插值出来的函数是这样的:
这就是龙格现象,主要原因是误差项里面的高阶导数导致的。什么是龙格现象(Rungephenomenon)?如何避免龙格现象?_MachineLearningwithTuring'sCat-CSDN博客_龙格现象www.baidu.comhttps://en.wikipedia.org/wiki/Runge%27s_phenomenonen.wikipedia.org
其实不单单是
,很多偶函数都有这样的性质:
比如
:
用于生成上图的代码:
# 龙格现象的产生
if __name__ == "__main__":
import numpy as np
import matplotlib.pyplot as plt
f = lambda x: 1/(1+25*x**2)
# 待插值的元素值
x_points = np.linspace(-1,1,11)
print(x_points)
y_points = list(map(f,x_points))
# 牛顿插值
x = np.linspace(-1,1)
y = list(map(lambda t: lagrange_interpolate(x_points,y_points,t),x))
# 画图
plt.figure("lagrange interpolation")
plt.scatter(x_points,y_points,color = "orange")
plt.plot(x,y)
plt.legend(["lagrange interpolation","scattered points"])
plt.show()
如何避免龙格现象,就需要用到分段插值了。
下一篇将讲解分段插值。
python插值程序_计算方法(2)——插值法(附Python程序)相关推荐
- python cnn 股市_股市分析——ATR指标(附python代码)
今天给大家介绍下ATR指标,废话不多说,直接上定义和公式,代码在公式后面(如果找不到文档的,可以从第一篇代码里下载) 纯理工男一枚,不喜勿喷,只是想编编代码玩玩,以后可能再用CNN.DNN什么的预测一 ...
- python 自动化微信小程序_干货 | 微信小程序自动化测试最佳实践(附 Python 源码)...
原标题:干货 | 微信小程序自动化测试最佳实践(附 Python 源码) 本文为霍格沃兹测试学院测试大咖公开课<微信小程序自动化测试>图文整理精华版. 随着微信小程序的功能和生态日益完善, ...
- 微信 小程序 python 渲染_干货 | 微信小程序自动化测试最佳实践(附 Python 源码)...
本文为霍格沃兹测试学院测试大咖公开课<微信小程序自动化测试>图文整理精华版,进阶学习文末加群! 随着微信小程序的功能和生态日益完善,很多公司的产品业务形态逐渐从 App 延升到微信小程序. ...
- 微信小程序python自动化测试_干货 | 微信小程序自动化测试最佳实践(附 Python 源码)...
本文为霍格沃兹测试学院测试大咖公开课<微信小程序自动化测试>图文整理精华版. 随着微信小程序的功能和生态日益完善,很多公司的产品业务形态逐渐从 App 延升到微信小程序.微信公众号等.小程 ...
- 微信小程序自动化测试最佳实践(附 Python 源码)
随着微信小程序的功能和生态日益完善,很多公司的产品业务形态逐渐从 App 延升到微信小程序.微信公众号等.小程序项目页面越来越多,业务逻辑也越来越复杂,全手工测试已无法满足快速增长的业务需求. 然而, ...
- 结束python服务器进程_服务器端后台持续执行python程序小demo
冰先生:python小脚本,爬天气预报并定时发邮件提醒(适合小情侣的甜蜜打开方式)zhuanlan.zhihu.com 上一篇写了做个可以发送邮件的小demo,这一篇呢我们来谈一谈如何让他更加方便的 ...
- python怎么做软件程序_如何打包和发布Python程序
如何打包和发布Python程序 在使用Python的过程中,我们经常需要做的一件事情就是通过pip来安装第三方的包.那么你是否也曾想过pip安装的包是怎么被打包并发布上去的呢?今天就来说一说Pytho ...
- python监控linux运行程序_如何在linux/tcl/python中监控正在打开或启动的应用程序?...
我正在尝试构建一个面板应用程序,类似于avant window navigator或UbuntuUnity.在 我的问题是,一旦我用预先定义好的应用程序构建了面板,当应用程序打开或启动时,如何向面板添 ...
- python恶搞图_搞几款由“Python”语言编写的“有趣、恶搞、好玩”的程序代码!...
分享一:"啥是佩奇?"让Python语言告诉你 用Python代码创作一副佩奇: # coding:utf-8 import turtle as t t.pensize(4) t. ...
- python嵌入到程序_在应用中嵌入Python:转
前面的章节讨论如何扩展Python,如何生成适合的C库等.不过还有另一种情况:通过将Python嵌入C/C++应用以扩展程序的功能.Python嵌入实现了一些使用Python更合适的功能.这可以有很多 ...
最新文章
- 69亿美元英伟达史上最大收购!这家基金又赢了
- 如何利用Gephi可视化浏览的网站关系
- NC:港大张彤团队-基于组学的耐药基因风险评估框架
- Dubbo 同步、异步调用的几种方式
- Python编程实现粒子群算法(PSO)详解
- python自动化测试视频教程_精品系列-悠悠Python自动化测试学习视频,资源教程下载...
- c语言编程抢30,抢三十-程序?
- python里的关键字有哪些_Python 中的关键字有哪些?
- 怎么调试EXC_BAD_ACCESS错误
- 监理公司的核心竞争力
- RHEL 6 配置yum源
- linux根文件系统的挂载过程详解
- 虚拟机python pip安装不了_给在Linux虚拟机里运行的FreeBSD 12安装pip Python包管理器...
- D*路径搜索算法原理解析及Python实现
- CAD虚线不显示怎么办?CAD虚线不显示解决办法
- jboot-admin学习
- Go简明语法汇总--入门
- Python Web前端实战案例——电商网站商品菜单导航栏
- 【Interview###】华为、中兴嵌入式(C)笔试题
- [vue] Vuex中四个map方法的使用 mapState mapGetters mapActions mapMutations
热门文章
- 算法导论-上课笔记7:贪心算法
- 一个程序员的江湖传奇!出江湖,闯江湖,退出江湖
- Linux每天一个命令:nc/ncat
- 计算两个数组之间重叠数字对的重叠个数(咋个办呢 zgbn)
- 2022年你应该知道的15 款 Python 编辑器/ IDE
- OSChina 周四乱弹 —— 举杯邀明月 对眼成三人
- 修改CentOS ll命令以 K 为单位显示文件大小
- 抖音卖货需要注意什么?传统电商如何入驻抖音带货
- Android 9.0 系统管控Wifi的启用和禁用功能实现
- ireport打印预览提示 The document has no pages的解决方法