Python数值计算(11)——拉格朗日插值

news2024/11/28 19:54:12

本篇介绍一下多项式插值中,拉格朗日法的原理及其实现。

1. 一点数学知识

先引用数学背景。如果给定N个点,然后要求一个多项式通过这N个点,最简单直接的方式是列出线性方程求解,N个点可以确定N个未知量,则所求的拟合多项式,其最高次幂就是(N-1)。

在很多教材上,关于该插值多项式的算法如下:

这个公式具有很强烈的轮换性,和线性方程组的克莱姆法则很相似(其实本质上就是一样的),所以实现起来也是非常简单。

2. 算法实现

def lagrangeIntp(x:ndarray,y:ndarray):
    '''
    返回基于x,y点对的拉格朗日插值多项式
    x,y必须有相同的长度
    '''
    n=len(x)
    ret=P([0])
    for k in range(n):
        roots=list(x)   
        s=roots[k]
        del roots[k]
        p=P.fromroots(roots)
        p=p*y[k]/p(s)
        ret+=p
    return ret

当然,上述函数并未对输入进行检测,例如x,y的尺寸应该是一致的等。

然后测试一下该拟合函数,仍旧以之前的函数f(x)=1-2x+3x^3为例:

a=P([1,-2,0,3])
print(a) # 1.0 - 2.0·x + 0.0·x² + 3.0·x³
x=np.arange(-2,2)
print(type(x))
y0=a(x)
b=lagrangeIntp(x,y0)
print(b) # 1.0 - 2.0·x + 0.0·x² + 3.0·x³
c=P.fit(x,y0,deg=3,domain=[-1,1])
print(c) # 1.0 - 2.0·x - (3.87289473e-15)·x² + 3.0·x³
X=np.linspace(-2,1,100)
plt.plot(X,a(X),'r')
plt.plot(X,b(X),'g--')
plt.plot(X,c(X),'b--')
plt.grid()
plt.show()

上面的例子中,先用原多项式a生成了四个点(-2,-19),(-1,0),(0,1),(1,2),然后将其作为拟合的输入,计算得到多项式b,并加入了通过fit拟合出来的多项式c,最后绘制三者的图形,如下:

可以看到三个图形是几乎完全重叠的,特别是使用拉格朗日法的插值多项式b,与原多项式a完全相同。

3. 现成工具

目前在numpy中暂未找到拉格朗日插值的模块,但是在Scipy中有对应的模块,那就是scipy.interpolate,用起来很简单,只要几行代码就可以搞定:

from scipy.interpolate import lagrange
import numpy as np
from numpy.polynomial import Polynomial as P
np.polynomial.set_default_printstyle('unicode')

x=np.array([-2,-1,0,1,2])
y=np.array([-19,0,1,2,21])
g=P(lagrange(x,y).coef[::-1])
print(g) 
# 1.0 - 2.0·x¹ + 1.1102230246251565e-16·x² + 3.0·x³

上面这个代码中,x,y是根据前面的函数生成的,注意,这里一共有五个点,所以最高次数有可能会达到4次,但是算法运行的结果基本上和前面的函数相同,只有x的2次幂有一个非常小的项,可以认为是浮点运算误差造成的。另外就是scipy.interpolate.lagrange返回的是一个numpy.poly1d的对象,这个在numpy已经是一个过时的类型了,语句P(lagrange(x,y).coef[::-1])是将其转化为了常用的Polynomial类型。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/1961404.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

下面关于枚举的描述正确的一项是?

A. 枚举中定义的每一个枚举项其类型都是String; B. 在Java中可以直接继承java.util.Enum类实现枚举类的定义; C. 利用枚举类中的values()方法可以取得全部的枚举项; D. 枚举中定义的构造方法只能够使用private权限声明; 答案选择…

springboot山东外事职业大学校园食堂点餐系统-计算机毕业设计源码10417

摘 要 近年来,随着国民收入的提高,各行业取得长足进步,也带动了互联网行业的快速发展,许多传统行业开始与互联网相结合,通过数字化转型打造新的发展生态。 本文针对山东外事大学校园食堂点餐系统的需求,基于…

Java内存区域与内存溢出异常详解

在Java编程中,理解Java虚拟机的内存布局及其管理机制对于开发高效、稳定的应用程序至关重要。Java虚拟机的内存主要分为几个运行时区域,这些区域各司其职,共同支撑起Java程序的运行。本文将详细探讨Java虚拟机的内存区域以及这些区域如何与内…

Yolov模型的使用及数据集准备(1)LabelImg的下载和使用

1、LabelImg下载: labelimg简单来说就是打标签用的软件,当需要使用自定义数据集进行模型训练时,往往需要使用该软件来打标签。 下载地址:GitHub - HumanSignal/labelImg 1.1下载之后对压缩包进行解压 2、打开电脑的anaconda pro…

MyBatis XML配置文件

目录 一、引入依赖 二、配置数据库的连接信息 三、实现持久层代码 3.1 添加mapper接口 3.2 添加UserInfoXMLMapper.xml 3.3 增删改查操作 3.3.1 增(insert) 3.3.2 删(delete) 3.3.3 改(update) 3.3.4 查(select) 本篇内容仍然衔接上篇内容,使用的代码及案…

8G 显存玩转书生大模型 Demo

创建可用环境 # 创建环境 conda create -n demo python3.10 -y # 激活环境 conda activate demo # 安装 torch conda install pytorch2.1.2 torchvision0.16.2 torchaudio2.1.2 pytorch-cuda12.1 -c pytorch -c nvidia -y # 安装其他依赖 pip install transformers4.38 pip in…

Moving Object Segmentation: All You Need Is SAM(and Flow) 论文详解

系列文章目录 文章目录 系列文章目录前言摘要1 引言2 相关工作3 SAM Preliminaries4 帧级分割Ⅰ:以流作为输入5 帧级分割Ⅱ:以流为提示6 序列级掩膜关联7 实验7.1 数据集7.2 评价指标7 .3 实施细节7.4 消融实验7.5 定量结果7 .定性可视化 8 结论致谢附录…

01 - 计算机组成原理与体系结构

文章目录 一,计算机系统硬件基本组成硬件软件 二,CPU的功能与组成功能组成运算器控制器 三,数据表示计算机的基本单位进制转换原码,反码,补码,移码数值表示范围浮点数表示 四,寻址五&#xff0c…

【Unity模型】古代亚洲建筑

在Unity Asset Store上,一款名为"Ancient Asian Buildings Pack"(古代亚洲建筑包)的3D模型资源包,为广大开发者和设计师提供了一个将古代亚洲建筑风格融入Unity项目的机会。本文将详细介绍这款资源包的特点、使用方式以…

如何选择合适的自动化测试工具!

选择合适的自动化测试工具是一个涉及多方面因素的决策过程。以下是一些关键步骤和考虑因素,帮助您做出明智的选择: 一、明确测试需求和目标 测试范围:确定需要自动化的测试类型(如单元测试、集成测试、UI测试等)和测试…

AI视频实战教程:DiffIR2VR-Zero-模糊视频8K高清修复技术

〔探索AI的无限可能,微信关注“AIGCmagic”公众号,让AIGC科技点亮生活〕 本文作者:AIGCmagic社区 猫先生 一、简 介 DiffIR2VR-Zero:一种创新的零样本视频恢复技术,该技术利用预训练的图像恢复模型,解决…

C++高性能通信:图形简述高性能中间件Iceoryx

文章目录 1. 概述2. 支持一个发布者多个订阅者2.2 Iceoryx为何不支持多个发布者发布到同一个主题 3. Iceoryx的架构和数据传输示意图3.1 发布者与订阅者的通信机制3.2 零拷贝共享内存通信机制 4. 使用事件驱动机制4.1 WaitSet机制4.2 Listener机制 5. 已知限制6. 参考 1. 概述 …

sci-hub下载不了的文献去哪里获取全文

我们在查找外文文献时经常会用到sci-hub,但sci-hub也有没有收录的文献,遇到这种情况我们可以用另一个途径来获取该文献。 例如这篇Wiley数据库中的文献:Unveiling Gating Behavior in Piezoionic Effect: toward Neuromimetic Tactile Sensin…

Linux服务管理(四)Apache服务

Apache服务 1、基于IP的虚拟主机2、基于IP端口的虚拟主机3、基于域名的虚拟主机4、prefork模式5、worker模式6、event模式7、细说驱动工作模式和MPM(多处理模块)工作模式 新旧域名都保留,因为旧域名已有一定的知名度和流量,直接下…

Cocos Creator2D游戏开发(8)-飞机大战(6)-炸机

碰撞 飞机与飞机碰撞 子弹与飞机碰撞 ① 设置碰撞矩阵 设置碰撞矩阵,就是设置谁跟谁碰撞(添加Enemy,PlayerBullet,Player) ②设置刚体和碰撞体 两个预制体设置(Enemy和PlayerBullet) 注意点: 1. 都在预制体节点上,不在图片上; 2.碰撞体Collider2D中的Editing悬着好之后可以调整…

C#-读取测序数据的ABI文件并绘制svg格式峰图-施工中

本地环境:win10,visual studio 2022 community 目录 前言问题描述解决思路实现效果 前言 本文是在已有的代码基础上进行的开发,前期已经实现: ABI文件的解析峰图的简单绘制svg绘图 对于1,主要用到之前重写的struct包…

大模型面经之bert和gpt区别

BERT和GPT是自然语言处理(NLP)领域中的两种重要预训练语言模型,它们在多个方面存在显著的区别。以下是对BERT和GPT区别的详细分析。 一、模型基础与架构 BERT: 全称:Bidirectional Encoder Representations from Trans…

系统移植(九)Linux内核移植(未整理)

文章目录 一、概念二、在linux内核源码的arch/arm/configs目录下生成FSMP1A板子对应的默认配置文件三、将自己编写的驱动通过图形化界面的方式编译到内核的镜像文件uImage中(一)拷贝myled.c和myled.h文件到linux内核源码的drivers/char目录下&#xff08…

第15周 15.1 Zookeeper简介安装及基础使用

1. Zookeeper介绍 1.1 介绍 1.2 应用场景简介 1.3 zookeeper工作原理 1.4 zookeeper特点

Canva收购Leonardo.ai,增强生成式AI技术能力

每周跟踪AI热点新闻动向和震撼发展 想要探索生成式人工智能的前沿进展吗?订阅我们的简报,深入解析最新的技术突破、实际应用案例和未来的趋势。与全球数同行一同,从行业内部的深度分析和实用指南中受益。不要错过这个机会,成为AI领…