数学建模——拟合与插值（插值部分）

最新推荐文章于 2022-07-27 00:28:33 发布

Blanche117

最新推荐文章于 2022-07-27 00:28:33 发布

阅读量3.4k

点赞数 9

分类专栏：数学建模文章标签： matlab 数值分析数值分析简明教程非线性拟合多元插值

本文链接：https://blog.csdn.net/weixin_45745854/article/details/107146815

版权

数学建模专栏收录该内容

12 篇文章 30 订阅

订阅专栏

插值与拟合

一、插值与拟合概念介绍

问题的引入
已经测得海洋的某深度的及其对应的水温：

在这里插入图片描述

如何根据这些已有的数据如何估计其他深度（比如600,700,800米处）的水温？我们很自然想到深度和水温之间是否存在某种函数关系。

在这里插入图片描述

函数的表达式可能无法给出，只能通过实验或者观察得到有限数量的数据点，那么如何通过数据点得到其他的点函数值？

在这里插入图片描述 2. 插值的概念：

在实际问题中，一个函数y=f(x)往往是通过实验观察到的，仅已知函数f(x)在某个区间[a,b]上一系列点的值
在这里插入图片描述
当需要这些节点

之间的某点x的函数值时，常用较为简单的，满足一定的条件的函数𝑃(𝑥)去代替真实的，难以得出的𝑓(𝑥)，插值法是一种常用的方法，其插值函数满足

拟合的概念
拟合也是有限个数据点，求近似函数。但是拟合只要求整体上逼近，而不要求一定满足①的条件，即不要求拟合得到的曲线一定过数据点，但是要求在某种意义上这些点的总偏差最小。
拟合和差值的比较

二、插值方法

Lagrange插值
① 插值问题描述
设函数 f(x)在区间[a,b]上有定义，在其n+1个互异点x_i∈[a,b]处的函数值 y_i= f(x_i) (i=0,1,…,n)为已知，求一个次数不超过n的代数多项式

使其满足

该问题成为插值问题。
其中：
xi：插值节点
P n (x)：拉格朗日插值多项式
多项式插值：P n (x)为多项式函数
分段插值：P n (x)为分段多项式函数
三角插值：P n (x)为三角函数
② 代数差值问题：
（1）满足Lagrange插值条件的多项式P_n (x)是否唯一？
有定理：对于n+1个互异的插值节点，Lagrange差值多项式存在且唯一。
证明：

系数行列式为范德蒙德行列式，由于插值节点互异，所以行列式不为零，所以有唯一的插值多项式存在。

（2）如何求P_n (x)，或者说构造出这个多项式？
定理：

满足：
定理：拉格朗日差值问题的解——

其中Lj (x)已由上面的定理给出。
例如当n=3时，上述的解即为：

（3） Lagrange插值编程
例：
求：

解：给出Matlab代码：

编写函数Lagrange

1.	function [waity,fnew]=lagrange(x0,y0,waitx)  
2.	n=length(x0);   %求出插值节点的数量,减1即是多项式的次数  
3.	m=length(waitx);    %求出待求点的数量  
4.	waity=zeros(1,m);  
5.	syms x temp  
6.	f=0;  
7.	for k=1:n  
8.	    temp=1;  
9.	    for j=1:n  
10.	        if j~=k  
11.	            temp=((x-x0(j))/(x0(k)-x0(j)))*temp;  
12.	        end  
13.	    end  
14.	    f=f+temp*y0(k);  
15.	end  
16.	fnew=collect(f,x);  
17.	for i=1:m  
18.	    waity(i)=subs(fnew,'x',waitx(i));  
19.	end  
20.	ezplot(fnew,[x0(1),x0(n)]); %绘制出插值函数的图像  
21.	hold on  
22.	plot(x0,y0,'*r');           %绘制出插值结点  
23.	hold on  
24.	plot(waitx,waity,'<b');     %绘制出待求点的值  
25.	legend('插值函数图像','插值节点','待求点');  
26.	title('Lagrange插值法插值函数图像');

命令行窗口输入代码：

1.	x0=[0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5];  
2.	y0=[0,0.4794,0.8415,0.9975,0.9093,0.5985,0.1411,-0.3508,-0.7568,-0.9775,-0.9589];  
3.	waitx=[0.6,1.2,2.8,3.7,5];  
4.	[waity,fnew]=lagrange(x0,y0,waitx)

绘制出的图像：
在这里插入图片描述
导出结果：

得到的插值函数为：

（4） Runge（龙格）现象简单介绍
Lagrange插值多项式的次数必须保持适度大小，过大或者过小都有可能导致较大的插值误差。
由于采用高次而产生较大误差的现象称为Runge (龙格)现象。
详细原因：https://www.zhihu.com/question/39329749
在这里插入图片描述
（5） Lagrange插值法的缺点：
无法增加采样点，采样点增加所有的系数需要重新计算。

分段低次插值
（1）将插值区间为[a, b]分为n个小区间[x_k,x_(k+1)], 在每一个小区间上使用低次插值多项式，其中，诸点（x_k,x_(k+1)）称为型值点，且y_k=f(x_k )。求过型值点（x_k,y_k）(k=0,1,2……n)的折线函数L_h (x)，满足以下条件：

（2）在Matlab中直接应用interp1(X, Y, x)或interp1(X, Y, x, ‘linear’) 命令即可。其中，X=[x0, x1, …, xn-1, xn]为节点向量, Y=[y0, y1, …, yn-1, yn]为对应的节点向量函数值，x为所求插值点。
（3）分段线性插值和Lagrange插值比较
下面给出Matlab中的比较绘图代码：

 1. function comparsionLH  
 2. X=[-5:5];  
 3. Y=cos(X)./sqrt(1+sin(X).^2);  
 4. x=[-5,-4.2,-3.5,-1,1.2,2.7,4.5];  
 5. [~,~]=lagrange(X,Y,x);  
 6. hold on  
 7. y=interp1(X,Y,x);  
 8. t=[-5:0.01:5];  
 9. f=cos(t)./sqrt(1+sin(t).^2);  
 10. plot(X,Y,'r*-',x,y,'b*',t,f,'k','markersize',10,'LineWidth',2.0)  
 11. title('Lagrange插值和分段低次插值比较图')  
 12. legend('插值函数图像','插值节点','待求点','分段低次插值曲线','低次样条插值插值节点','真实曲线')

代码生成图片
由上图可知分段低次插值的效果比拉格朗日插值的效果要好，更加逼近真实情况：

Newton插值公式
（1）差商的定义（略）
（2）公式：
由数学归纳法可得n阶差商公式：

（3） Newton插值公式

（4）优点：差商可以由“同列威尔法”填表得到，在增加采样点时易于操作，不用重新计算所有的系数。
同列威尔法：

（5） Newton插值公式程序实现：（仍以上述Lagrange插值法的表格数据为例）
p.s.导入插值节点的命令行代码略去，与上述的拉格朗日插值法相同，可以直接保存为mat文件，便于后续代码操作。

4.	function [y,fnew,d]=newton_interpolation(x0,y0,waitx)  
5.	n=length(x0);  
6.	m=length(waitx);  
7.	d=zeros(n-1,n-1);  
8.	y=zeros(m,1);  
9.	syms x  
10.	%模拟填表过程  
11.	for k=1:n-1               
12.	    d(k,1)=(y0(k+1)-y0(k))/(x0(k+1)-x0(k)); %求一阶差商，即准备好初始值  
13.	end  
14.	%模拟从第二列开始填表  
15.	for k=2:n-1  
16.	    for j=2:k  
17.	        d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x0(k+1)-x0(k+1-j));  
18.	    end  
19.	end  
20.	%求出系数向量a  n维列向量  
21.	a=[y0(1);diag(d)];  
22.	  
23.	%求出函数表达式  
24.	f=y0(1);  
25.	temp=1;  
26.	for i=2:n  
27.	    for j=1:i-1  
28.	        temp=temp*(x-x0(j));  
29.	    end  
30.	    f=f+a(i)*temp;  
31.	    temp=1;  
32.	end  
33.	fnew=collect(f,x);  
34.	for k=1:m  
35.	    y(k)=subs(fnew,'x',waitx(k));  
36.	end  
37.	plot(x0,y0,'-.b',waitx,y,'r*');  
38.	hold on  
39.	ezplot(f,[x0(1),x0(n)]);  
40.	grid on  
41.	title('Newton插值法函数图像');  
42.	legend('插值节点连线','待求节点','Newton插值法插值曲线');

在这里插入图片描述

三、Matlab插值工具箱

一维插值函数
Matlab中有现成的一维插值函数interp1，语法如下：

y=interp1(x0,y0,x,‘method’)
其中：method指定的插值的方法默认为线性插值。其值可为：
‘nearest’------最近项插值
‘linear’------线性插值
‘spline’------立方样条插值
‘cubic’------立方差值
所有的插值方法都要求x0单调
当x0等距时可以使用快速插值法，使用快速插值法的格式为’*nearest ’ , ’ *linear’,’*spline’,’*cubic’

三次样条插值

在Matlab中数据点成为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结点(not-a-knot)条件。这个条件强迫第一个和第二个三次多项式的三阶导数相等。对于最后一个和倒数第二个三次多项式也做同样的处理。
Matlab中的三次样条插值有以下的函数：
y=interp1(x0,y0,waitx,‘spline’)；
y=spline(x0,y0,waitx);
pp=csape(x0,y0,conds);
pp=csape(x0,y0,conds,valconds);y=fnval(pp,waitx);

其中x0,y0是已知数据点;waitx是待求的插值点，y是插值点的返回值
对于三次样条插值，提倡使用函数caspe。csape的返回值是pp形式，要求插值点的函数值必须调用函数fnval
上述函数的参数中的’conds’即边界条件。
如果没有写的话即采用默认的边界条件：Lagrange边界条件
以下介绍几种边界条件：

'not-a-knot’非扭结点条件
'periodic’周期条件
'second’边界为二阶导数，二阶导数的值在valconds参数中给出，若忽略valconds参数则二阶导数的默认值是[0, 0]

四、例题

待加工零件的外形根据工艺要求由一组数据(x,y)给出（在平面情况下），用机床加工时每一刀只能沿着x方向和y方向走很小的一步，这就需要从已知数据得到加工所要求的步长很小的(x,y)坐标。

下表给出的x,y数据位于机翼断面的下轮廓线上，假设需要得到x的坐标每改变0.1时的y的坐标。试完成加工所需要的数据，画出曲线，并求出x=0时的曲线斜率和13≤x≤15范围内的y的最小值。要求使用分段线性和三次样条两种插值方法进行计算。
在这里插入图片描述
解：首先给出matlab实现的代码：

1.	function [xmin,ymin,dy_dx0,data]=myMachineFinish  
2.	x0=[0   3   5   7   9   11  12  13  14  15];  
3.	y0=[0   1.2  1.7    2.0  2.1    2.0  1.8    1.2  1.0    1.6];  
4.	x=0:0.1:15;                     %要求的步长为0.1  
5.	y1=interp1(x0,y0,x);            %采用默认的分段线性插值  
6.	y2=interp1(x0,y0,x,'spline');   %采用三次样条插值的spline方法  
7.	pp1=csape(x0,y0,'second');      %采用三次样条插值的second方法，二阶导数的默认值为[0,0]；  
8.	y3=fnval(pp1,x);  
9.	pp2=csape(x0,y0);               %采用默认的拉格朗日边界条件  
10.	y4=fnval(pp2,x);  
11.	data=[x',y1',y2',y3',y4'];  
12.	subplot(1,3,1)  
13.	plot(x0,y0,'+',x,y1)            %画出给出的数据点以及采用分段线性插值的结果  
14.	title('Piecewise linear')  
15.	subplot(1,3,2)  
16.	plot(x0,y0,'+',x,y2)            %画出给出的数据点以及采用三次样条插值的结果  
17.	title('spline1');  
18.	subplot(1,3,3)  
19.	plot(x0,y0,'+',x,y4)            %画出给出的数据点以及采用三次样条插值csape函数的结果  
20.	title('spline2');  
21.	dx=diff(x);  
22.	dy=diff(y4);  
23.	dy_dx=dy./dx;                   %求出要求的步长点的导数  
24.	dy_dx0=dy_dx(1);                %求出x=0处的斜率  
25.	ytemp=y4(131:151);              %取出x∈[13,15]之间的所有y的值  
26.	ymin=min(ytemp);                %求出最小值  
27.	index=find(y4==ymin);           %找到最小值的索引  
28.	xmin=x(index);                  %找到对应索引的x的值

命令行窗口代码：

[xmin,ymin,dy_dx0,data]=myMachineFinish

输出结果：（data部分数据过多仅仅展示一部分）

1.	xmin =  
2.	   13.8000  
3.	ymin =  
4.	    0.9851    
5.	dy_dx0 =  
6.	    0.4972    
7.	data =  
8.	         0         0         0         0         0  
9.	    0.1000    0.0400    0.0499    0.0441    0.0497  
10.	    0.2000    0.0800    0.0990    0.0881    0.0987  
11.	    0.3000    0.1200    0.1474    0.1321    0.1470  
12.	    0.4000    0.1600    0.1951    0.1760    0.1946  
13.	    0.5000    0.2000    0.2421    0.2198    0.2415  
14.	    0.6000    0.2400    0.2884    0.2635    0.2878  
15.	    0.7000    0.2800    0.3340    0.3070    0.3333  
16.	    0.8000    0.3200    0.3788    0.3503    0.3781  
17.	    0.9000    0.3600    0.4230    0.3934    0.4223  
18.	    1.0000    0.4000    0.4665    0.4362    0.4658  
19.	    1.1000    0.4400    0.5094    0.4788    0.5086

在这里插入图片描述

Blanche117

关注

9
点赞
踩
54

收藏

觉得还不错? 一键收藏
打赏
0
评论
数学建模——拟合与插值（插值部分）

插值与拟合一、插值与拟合概念介绍问题的引入已经测得海洋的某深度的及其对应的水温：如何根据这些已有的数据如何估计其他深度（比如600,700,800米处）的水温？我们很自然想到深度和水温之间是否存在某种函数关系。函数的表达式可能无法给出，只能通过实验或者观察得到有限数量的数据点，那么如何通过数据点得到其他的点函数值？ 2. 插值的概念：在实际问题中，一个函数y=f(x)往往是通过实验观察到的，仅已知函数f(x)在某个区间[a,b]上一系列点的值当需要这些节点之间的某点x的函数
复制链接

扫一扫