数据插值之朗格朗日插值(一)

目录

一、引言

二、代码实现 

2.1 Lagrange插值求插值多项式:

代码解析:

1.vpa解释

 2.ploy(x)解释:

3.conv()解释

4.poly2sym()解释

2.2 Lagrange插值求新样本值和误差估计:

代码解析:

1.errorbar(x, y, R, '.g')

2.plot(X, Y, 'or')

三、Lagrange在Matlab中的使用模板


一、引言

数据插值是指通过有限个原始数据点构造出一个解析表达式,从而可以计算数据点之间的函数值。在Matlab中数据插值的方法主要有拉格朗日插值、牛顿差值、三次样条插值、埃尔米特插值、一维数据插值、二维数据插值等等。笔者旨在写一个使用Matlab进行数值插值的专题,通过介绍插值原理,代码实现,形成模板,方便日后使用时直接代入参数直接调用。

本节介绍拉格朗日插值(Lagrange interpolation)。Lagrange插值是一种用于逼近通过给定数据点的函数的方法。它涉及构造一个度数为 ( n-1 ) 或更低的多项式,该多项式通过 ( n ) 个给定点。该多项式是使用Lagrange基多项式构造的,这是一组多项式,当在给定点之一进行评估时,会得到1,并且当在任何其他给定点进行评估时会得到0。

Lagrange插值多项式可表示为:

L(x) = \sum_{i=0}^{n} y_i l_i(x)

Lagrange插值基函数可以表示为:

l_i(x) = \frac{(x-x_0)(x-x_1)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i - x_0)(x_i - x_1)\cdots(x_i - x_{i-1})(x_i - x_{i+1})\cdots(x_i - x_n)}

二、代码实现 

2.1 Lagrange插值求插值多项式:

function[C, L, L1, l] = lagranl(X, Y)
% 定义拉格朗日插值多项式和基函数
%
% 输入:
%   X - 插值节点的横坐标向量
%   Y - 插值节点的纵坐标向量
%
% 输出:
%   C - 拉格朗日插值多项式的系数向量
%   L - 拉格朗日插值多项式的符号表达式
%   L1- 拉格朗日基函数对应的系数矩阵
%   l - 拉格朗日基函数对应的符号表达式

m = length(X); % 获取插值节点的个数

L1 = zeros(m, m); % 初始化存储拉格朗日基函数系数的矩阵
l = sym(zeros(1, m)); % 初始化单行符号值存储基函数表示式

for k = 1:m % 遍历每一个插值节点
    V = 1; % 初始化基函数的系数向量为1
    for i = 1 : m % 遍历所有节点用于计算
        if k ~= i %排除掉自身的情况
            V = conv(V, poly(X(i))) / (X(k) - X(i)); % 更新基函数的系数向量
        end
    end
    L1(k, :) = V; % 存储第k个基函数的系数
    l(k) = poly2sym(V) * Y(k); % 将基函数乘以对应的Y值,转换为符号表达式并存储
end

C = sum(Y .* L1, 1); % 计算拉格朗日插值多项式的系数向量
L = sum(l); % 计算拉格朗日插值多项式的符号表达式

 下面进行测试:

X = [-2.2, -1.0, 0.01, 1.0, 2.0, 3.3, 2.2];
Y = [17.1, 7.3, 1.1, 2.0, 17.1, 23.1, 19.3];
[C, L, L1, l] = lagranl(X, Y);
L = vpa(L, 3)

代码解析:

1.vpa解释

vpa 是 MATLAB 中的一个函数,表示 “Variable Precision Arithmetic”(可变精度计算)。在这个特定的调用中,vpa(L, 3) 将 L 的符号表达式以3位小数的精度表示。具体解释:

  • L: 是从 lagranl 函数中返回的拉格朗日插值多项式的符号表达式。
  • vpa(L, 3): 会将这个符号表达式的系数和结果值保留到小数点后3位。

假如:

L = sym('0.123456789*x^2 + 0.987654321*x + 3.141592653');

使用 vpa 后:

L = vpa(L, 3)
% 将变为
L = sym('0.123*x^2 + 0.988*x + 3.14');
 2.ploy(x)解释:

在MATLAB中,ploy(x)函数并不是直接用来创建多项式的,容易与用于多项式插值的ployfit()或构造多项式系数的polyval()函数混淆。实际上,ploy(x)函数是用来从根求多项式的系数的。也就是说,如果你有一系列的复数或实数根,poly() 可以帮助你找到一个多项式,其在这些点上的值为零。

 基本用法:

p = poly(r)
  • r 是一个向量,包含了多项式的根。
  • p 是返回的结果,是一个向量,表示对应于根 r 的多项式的系数,按降幂排列。

举个例子:

r = [1, 2];
p = poly(r);
disp(p);  % 输出多项式系数

 这段代码会输出 [1 -3 2],对应于多项式x^2-3x+2,验证了1和2确实是这个多项式的根。

3.conv()解释

conv是MATLAB 中用来执行多项式卷积(Convolution)的函数。在多项式运算中,conv 可以用来计算两个多项式的乘积。

 基本用法

C = conv(A, B)

其中 A 和 B 是两个多项式的系数向量,返回向量 C 表示这两个多项式的乘积的系数。

示例:

假设我们有两个多项式:

P(x) = 1 + 2x + 3x^2

Q(x) = 4 + 5x

用系数向量表示:

P = [3, 2, 1]; % 对应:1 + 2x + 3x^2
Q = [5, 4]; % 对应: 5x + 4

我们使用conv来计算两者的乘积:

P = [3, 2, 1];  % 3x^2 + 2x + 1
Q = [5, 4];     % 5x + 4

C = conv(P, Q)

% 输出 C
% C = [15 22 13 4] 对应最终方程式:15x^3 + 22x^2 + 13x + 4

4.poly2sym()解释

poly2sym()是把多项式系数转换为符号多项式。

V = [1, 2, 3]\rightarrow x^2+3x+2

2.2 Lagrange插值求新样本值和误差估计:

%% 拉格朗日插值及误差估计
% 此函数使用拉格朗日插值公式来计算插值值,并估计误差
% 输入参数:
% X - 已知数据点的x坐标 (向量)
% Y - 已知数据点的y坐标 (向量)
% x - 需要插值的x坐标 (向量)
% M - 最大连续(n+1)阶导数的上界 (标量)
% 输出参数:
% y - 插值后的y值 (向量)
% R - 误差估计 (向量)

function [y, R] = lagranzi(X, Y, x, M)
    n = length(X);      % 已知数据点的数量
    m = length(x);      % 需要插值的数据点数量

    for i = 1:m
        z = x(i);       % 当前需要插值的x坐标
        s = 0.0;        % 插值和初始化为0

        for k = 1:n
            p = 1.0;    % 插值多项式L_k(z)初始化为1
            q1 = 1.0;   % 误差估计中的分子初始化为1
            c1 = 1.0;   % 误差估计中的分母初始化为1

            for j = 1:n
                if j ~= k
                    % 计算拉格朗日基函数L_k(z)
                    p = p * (z - X(j)) / (X(k) - X(j));
                end
                % 计算误差估计的分子部分
                q1 = abs(q1 * (z - X(j)));
                % 计算误差估计的分母部分
                c1 = c1 * j;
            end
            
            % 加权求和得到插值值
            s = p * Y(k) + s;
        end
        
        y(i) = s;         % 存储插值结果
        R(i) = M * q1 / c1;  % 计算并存储误差估计
    end
end

下面进行测试:

clc
clear
X = [-2.2, -1.0, 0.01, 1.0, 2.0, 3.3, 2.2];
Y = [17.1, 7.3, 1.1, 2.0, 17.1, 23.1, 19.3];
x = linspace(-3, 4, 50);
M = 1;
[y, R] = lagranzi(X, Y, x, M);
errorbar(x, y, R, '.g')
hold on
plot(X, Y, 'or')
x = 2.8;
[y, R] = lagranzi(X, Y, x, M);
x = -3:0.01:4;
L = 0.21*x.^6 - 0.87*x.^5 - 1.21*x.^4 + 5.9*x.^3 + 4.48*x.^2 - 7.68*x + 1.18;
plot(x, L)
legend('误差', '样本点', '拉格朗日多项式函数曲线')
print(gcf, '-r600', '-djpeg', '图2-1.jpg')

代码解析:

1.errorbar(x, y, R, '.g')

在MATLAB中,errorbar() 函数用于绘制带有误差条的图形,它能够直观地展示数据点的不确定度或误差范围。函数调用格式中的各个参数含义如下:

errorbar(x, y, R, '.g')

这里各参数的含义是:

  • x:数据点,在x轴上的值。
  • y:数据点,与x对应,在y轴上的值。
  • R:这可以是一个向量或一个矩阵,定义了误差条的长度。如果是向量,那么这个向量提供了每个数据点y方向上的误差估计;如果是矩阵(通常是2列),第一列是负方向的误差,第二列是正方向的误差,用于分别表示y值的下限和上限。
  • '.':这是一个标记符号参数,指定了数据点的样式,在这个例子中使用的是点。
  • 'g':颜色代码,指定了数据点和误差条的颜色,在这里是绿色。
2.plot(X, Y, 'or')
  • 'o' 表示数据点将以圆圈形式标记。如果只用 'o',MATLAB 会在每个 (X, Y) 数据对的位置绘制一个空心的圆圈。
  • 'r' 表示红色(red)。结合前面的 'o',这意味着每个数据点将以红色填充的圆圈来表示。

三、Lagrange插值使用模板

如果目前你拥有X,Y样本点,仅想求出一个新的x所对应的y值,可以直接在下述代码中进行修改“示意用法”中的值,程序放入Matlab中可以直接执行:

clc
clear
% 示例用法
x = [1, 2, 3, 4];
y = [1, 4, 9, 16];
xi = 2.5;
L_value = lagrange_interpolation(x, y, xi);
disp(['在 xi = ' num2str(xi) ' 处的插值值是 ' num2str(L_value)]);

function L = lagrange_interpolation(x, y, xi)
    % 拉格朗日插值函数
    % x 和 y 是已知的数据点
    % xi 是希望插值的点
    % 返回的 L 是在 xi 处的插值结果
    
    % 确保 x 和 y 长度一致
    if length(x) ~= length(y)
        error('x 和 y 必须具有相同的长度');
    end
    
    % 初始化结果 L
    n = length(x);
    L = 0;
    
    % 计算拉格朗日多项式
    for i = 1:n
        % 初始化基多项式
        Li = 1;
        for j = 1:n
            if i ~= j
                Li = Li * (xi - x(j)) / (x(i) - x(j));
            end
        end
        % 累加到结果
        L = L + Li * y(i);
    end
end

参考资料:《Matlab编程与汽车仿真应用》 ——崔胜民

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

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

容器化:ES和Kibana

1 缘起 最近在学习使用ES, 为了找一个功能强大的可视化工具,之前使用了ES-Head,可以满足学习需求。 闲暇时间又折腾了另一个工具Kibana, 分享如下。 Kibana优点: 用户友好性:Kibana提供直观易用的用户界面…

项目思考-编辑器

1、文本生成编辑器 2、图片合成编辑器(未完待续) 3、文字和图像版本的技术要点,区别(未完待续) 4、编辑器的人员配置考虑,技术难点分析(未完待续) 1、文本生成编辑器

【Python】 掌握 Flask 请求数据获取的艺术

基本原理 在Web开发中,Flask是一个用Python编写的轻量级Web应用框架。它被广泛用于快速开发简单的Web应用。当用户通过浏览器或其他客户端向服务器发送请求时,Flask需要能够接收和解析这些请求中的数据。这些数据可以是GET请求的查询字符串、POST请求的…

收集 VSCode 常用快捷键

快速复制行 Shift Alt ↑/↓ 都是往下复制行,区别是:按↓复制时光标会跟着向下移动,按↑复制时光标不移动。 向上/向下移动一行 Alt ↑/↓ 删除整行 Ctrl Shift KCtrl x 剪切快捷键在VSCode 可以直接删除一行 垂直编辑 Ctrl…

用于时间序列概率预测的蒙特卡洛模拟

大家好,蒙特卡洛模拟是一种广泛应用于各个领域的计算技术,它通过从概率分布中随机抽取大量样本,并对结果进行统计分析,从而模拟复杂系统的行为。这种技术具有很强的适用性,在金融建模、工程设计、物理模拟、运筹优化以…

HACL-Net:基于MRI的胎盘植入谱诊断的分层注意力和对比学习网络

文章目录 HACL-Net: Hierarchical Attention and Contrastive Learning Network for MRI-Based Placenta Accreta Spectrum Diagnosis摘要方法实验结果 HACL-Net: Hierarchical Attention and Contrastive Learning Network for MRI-Based Placenta Accreta Spectrum Diagnosis…

Linux驱动设备导论(1)

最近本人在学习Linux驱动,本系列教程是本人在一边学习,一边总结的系列教程,希望能够给很多刚学驱动小伙伴一些总结。 1.Linux设备分类 驱动针对的对象是存储器和外设,不是针对CPU,可以分为以下三大类: 1.…

01-Linux【准备篇】

一、学Linux的作用? 1.Linux下开发(部署)软件项目 2.Linux运维 二、Linux的强与弱 1.薄弱 个人桌面领域的应用 此领域是传统Linux应用薄弱的环节,近些年随着Ubuntu、fedora等优秀桌面环境的兴起,Linux在个人桌面领域的占有率在慢慢提高…

[国产大模型简单使用介绍] 开源与免费API

个人博客:Sekyoro的博客小屋 个人网站:Proanimer的个人网站 随着大模型技术蓬勃发展和开源社区越来越活跃,国内的大模型也如雨后春笋一般.这时,一些就会问了,有了llama3,Mistral还有Gemma等等,国外大厂接连发力,一些开源社区也会有一些不错的模型,国内怎么比?对一个人使用,oll…

Debezium+Kafka:Oracle 11g 数据实时同步至 DolphinDB 解决方案

随着越来越多用户使用 DolphinDB,各式各样的应用场景对 DolphinDB 的数据接入提出了不同的要求。部分用户需要将 Oracle 11g 的数据实时同步到 DolphinDB 中来,以满足在 DolphinDB 中实时使用数据的需求。本篇教程将介绍使用 Debezium 来实时捕获和发布 …

03_前端三大件CSS

文章目录 CSS用于页面元素美化1.CSS引入1.1style方式1.2写入head中,通过写style然后进行标签选择器加载样式1.3外部样式表 2.CSS样式选择器2.1 元素选择器2.2 id选择器2.3 class选择器 3.CSS布局相关3.1 CSS浮动背景:先设计一些盒子因此,引出…

【Go专家编程——内存管理——垃圾回收】

垃圾回收 所谓的垃圾就上不在需要的内存块,垃圾如果不清理,这些内存块就没有办法再次被分配使用。在不支持垃圾回收的编程语言中,这些垃圾内存就上泄露的内存。 1. 垃圾回收算法 常见的垃圾回收算法有3种 引用计数:对每个对象…

Vue学习笔记3——事件处理

事件处理 1、事件处理器(1)内联事件处理器(2)方法事件处理器 2、事件参数3、事件修饰符 1、事件处理器 我们可以使用v-on 指令(简写为)来监听DOM事件,并在事件触发时执行对应的JavaScript。 用法: v-on:click"me…

牛客NC334 字典序第K小【困难 10叉树 Java/Go/PHP/C++】,力扣 440. 字典序的第K小数字

题目 题目链接: https://www.nowcoder.com/practice/670c2bda374241d7ae06ade60de33e8b https://leetcode.cn/problems/k-th-smallest-in-lexicographical-order/description/ 本答案核心 10叉树, 数学规律Java代码 import java.util.*;public class Solution {…

出题123

题目时限空间说明 无特殊均默认 1 s , 256 M B 1s,256MB 1s,256MB Problem a 最大化 在最大化目标值的基础上选择的操作越多越好,且输出操作应当按照顺序执行,即你的输出顺序就是你的执行顺序,当有多个执行顺序可以最大化目标值时&#xff0…

49 序列化和反序列化

本章重点 理解应用层的作用,初识http协议 理解传输层的作用,深入理解tcp的各项特性和机制 对整个tcp/ip协议有系统的理解 对tcp/ip协议体系下的其他重要协议和技术有一定的了解 学会使用一些网络问题的工具和方法 目录 1.应用层 2.协议概念 3. 网络计…

网络爬虫原理及其应用

你是否想知道Google 和 Bing 等搜索引擎如何收集搜索结果中显示的所有数据。这是因为搜索引擎对其档案中的所有页面建立索引,以便它们可以根据查询返回最相关的结果。网络爬虫使搜索引擎能够处理这个过程。 本文重点介绍了网络爬虫的重要方面、网络爬虫为何重要、其…

Docker学习(3):镜像使用

当运行容器时,使用的镜像如果在本地中不存在,docker 就会自动从 docker 镜像仓库中下载,默认是从 Docker Hub 公共镜像源下载。 一、列出镜像列表 可以使用 docker images 来列出本地主机上的镜像。 各个选项说明: REPOSITORY&am…

vue源码2

vue之mustache库的机理其实是将模板字符串转化为tokens 然后再将 tokens 转化为 dom字符串&#xff0c;如下图 对于一般的将模板字符串转化为dom字符串&#xff0c;这样不能实现复杂的功能 let data {name:小王,age:18 } let templateStr <h1>我叫{{name}},我今年{{ag…

MySQl创建数据库与管理表

创建数据库与管理表 基础知识 完整的数据存储过程 同时&#xff0c;数据库系统层次 数据库服务器 -》 数据库 -》 数据表 -》 行与列 数据库命名规则&#xff1a; 库名、表名不得超过30字符&#xff1b;变量名&#xff08;字段&#xff09;不超过29字符 只能包含A-Z、a-z、…