当前位置:网站首页>Lagrange插值公式matlab实现
Lagrange插值公式matlab实现
2022-08-09 15:57:00 【泡泡怡】
一、公式推导原理
N次插值基函数:
满足插值多项式
形如此公式的插值多项式称为Lagrang插值多项式。
由的定义可知
若引入计号,再求导
因此
二 、符号说明
输入:
xi:已知数据点的横坐标
k:函数lk(x)的下标k
xx:待插值点的横坐标
输出:lk_x,即函数lk(x)在xx坐标点的纵坐标
代码如下:
function li_x = LagrangeFactor( xi, i, xx )
w = 1;
n = length( xi );
syms x;
for j = 1 : n
w = w * ( x - xi(j) );
end
dw = diff( w );
dwf = matlabFunction( dw );
dwi = dwf( xi(i) );
lx = ( w / ( x - xi(i) ) ) / dwi;
f = matlabFunction( lx );
lk_x = f( xx );
end
三、一次插值
1. 自变量函数的准备工作:
(xi,yi):是已知的数据点坐标
代码:xi = [ 0, 1 ];
yi = sin( xi );
n = length( xi );
y = 0;
x = [ xi(1) - 1 : 0.1 : xi(2) + 1 ];
2.根据lagrange插值多项式计算x坐标点处的函数值(纵坐标)
代码:
for k = 1 : n
lkx = LagrangeFactor( xi, k, x );
y = y + yi(k) * lkx;
end
3.绘图
代码:
figure;
plot( xi, yi, 'b.', 'markersize', 30 )
hold on
plot( x, sin(x), 'k--', 'LineWidth', 1.5 )
plot( x, y, 'r-', 'LineWidth', 2 )
legend( '插值点', '原曲线', '插值多项式曲线' );
axis( [ -1, 2, -1, 3 ] )
结果如图:
四、抛物插值
1.自变量函数的准备工作:
xi = [ 0, pi/2, pi ];
yi = [ 0, 1, 0 ];
n = length( xi );
y = 0;
x = [ xi(1) - 1 : 0.1 : xi(n) + 1 ];
2.根据lagrange插值多项式计算x坐标点处的函数值
for k = 1 : n
lkx = LagrangeFactor( xi, k, x );
y = y + yi(i) * lkx;
end
3.绘图
plot( xi, yi, 'b.', 'markersize', 30 )
hold on
plot( x, sin(x), 'k--', 'LineWidth', 1.5 )
plot( x, y, 'r-', 'LineWidth', 2 )
legend( '插值点', '原曲线', '插值多项式曲线' );
axis( [ xi(1) - 1, xi(n) + 1, -1, 2 ] )
汇总代码:
clear all
clc
%% 一次插值
%(xi,yi):
xi = [ 0, 1 ];
yi = sin( xi );
n = length( xi );
y = 0;
x = [ xi(1) - 1 : 0.1 : xi(2) + 1 ];
for k = 1 : n
lkx = LagrangeFactor( xi, k, x );
y = y + yi(k) * lkx;
end
%y
figure;
plot( xi, yi, 'b.', 'markersize', 30 )
hold on
plot( x, sin(x), 'k--', 'LineWidth', 1.5 )
plot( x, y, 'r-', 'LineWidth', 2 )
legend( '插值点', '原曲线', '插值多项式曲线' );
axis( [ -1, 2, -1, 3 ] )
%% 抛物插值
clear all
clc
xi = [ 0, pi/2, pi ];
yi = [ 0, 1, 0 ];
n = length( xi );
y = 0;
x = [ xi(1) - 1 : 0.1 : xi(n) + 1 ];
for k = 1 : n
lkx = LagrangeFactor( xi, k, x );
y = y + yi(k) * lkx;
end
%y
figure;
plot( xi, yi, 'b.', 'markersize', 30 )
hold on
plot( x, sin(x), 'k--', 'LineWidth', 1.5 )
plot( x, y, 'r-', 'LineWidth', 2 )
legend( '插值点', '原曲线', '插值多项式曲线' );
axis( [ xi(1) - 1, xi(n) + 1, -1, 2 ] )
function lk_x = LagrangeFactor( xi, k, xx )
w = 1;
n = length( xi );
syms x;
for j = 1 : n
w = w * ( x - xi(j) );
end
dw = diff( w );
dwf = matlabFunction( dw );
dwi = dwf( xi(k) );
lx = ( w / ( x - xi(k) ) ) / dwi;
f = matlabFunction( lx );
lk_x = f( xx );
end
边栏推荐
- Became CTO, was killed by my boss in 6 months, I lost 10 million
- 网络——数据交换方式
- 硬件开发的发展前景
- 网络——IPV4地址(二)
- Leetcode 算法面试冲刺 热题 HOT 100 刷题(406 416 437 438 448)(六十九)
- 【燃】是时候展现真正的实力了!一文看懂2022华为开发者大赛技术亮点
- 2022年中国第三方证券APP创新专题分析
- 二分法
- shopee引流方式有哪些,商家如何为自己店铺做引流?
- 1.1、VIFB: A Visible and Infrared Image Fusion Benchmark(一个可见光与红外图像融合Benchmark)文章阅读
猜你喜欢
随机推荐
网络——涉及的相关协议和设备汇总
B45 - 基于STM32单片机的家庭防火防盗系统的设计
CryoEM粒子(Particle)类型预测的数据集构建
QuickSort(快速排序)&&MergeSort(归并排序)的效率比较[搭配LeetCode例题]
MySQL索引的B+树到底有多高?
The Chinese Academy of Sciences slaps Google in the face: ordinary computers catch up with quantum superiority, and can solve calculations that would have taken 10,000 years in a few hours...
2022年8月9日:用C#生成.NET应用程序--使用 Visual Studio Code 调试器,以交互方式调试 .NET 应用(不会,失败)
PADS generates bitmap
Now, how to choose a stage rental LED display?
总结了 110+ 公开专业数据集
价值10亿美元 美国向乌克兰提供单次最大规模安全援助
网络——虚拟专用网和地址转换NAT
微信开发者工具程序开发好后,不报错,但是黑屏「建议收藏」
Knowledge Bits - How to Write a Project Summary
Arrow parquet 之 String Reader
SQL trill interview: send you a universal template, to?(key, each user to log on to the maximum number of consecutive monthly)
【教程3】疯壳·ARM功能手机-整板资源介绍
计组——大端方式和小端方式相关题目
安装MySQL的最后一步第四个红叉怎么解决
ceph部署