TestDoc

这是一篇测试文章，选用的内容为数学建模插值的内容：

插值

[toc]

一、插值定义

数据处理问题：在平面上给定一组离散点列，要求一条曲线把这些点依次连接起来，称之为插值

注意！此部分出现了大量不同的有关矩阵的函数运算方法，不同的运算方法之间对行标和列标的对应不同，最好在使用之前用size()方法看清相应的大小和对应关系，再进行操作！

二、常见插值方法

已知 $n + 1$ 个点 $(x_{i}, y_{i}) (i = 0, 1, \dots, n)$ ，下面求各种插值函数

（一）分段线性插值

将两个相邻的节点用直线连接起来，如此形成的一条折线就是分段线性插值函数，记作$In(x) $，它满足$ I_n(x_i)=y_i $，且$ I_n(x) $在每一个小区间$ [x_i,x{i+1}] $上是线性函数$ (i=0,1,…,n)$

$In(x) $可以表示为$ I_n(x)\sum\limits^n{i=0}y_il_i(x)$，其中

l_{i} (x) = {\begin{cases} \frac{x - x_{i - 1}}{x_{i} - x_{i - 1}}, x \in [x_{i - 1}, x_{i}] (i \neq 0) \\ \frac{x - x_{i + 1}}{x_{i} - x_{i + 1}}, x \in [x_{i}, x_{i + 1}] (i \neq n) \\ 0, O t h e r s \end{cases}

$I_{n} (x)$ 有较好的收敛性，即对于 $x \in [a, b]$ ，有

lim_{n \to \infty} I_{n} (x) = f (x)

用 $I_{n} (x)$ 计算x点的插值时，只用到x左右的两个节点，计算量与节点个数n无关。但n越大，分段越多，插值误差越小。实际上用函数表做插值计算时，分段线性插值就足够了，如：数学、物理重用的特殊函数表，数理统计中用的概率分布表等等

（二）拉格朗日插值多项式

拉格朗日插值的基函数为

\begin{aligned} (1) & l_{i} (x) & = \frac{(x - x_{0}) \dots (x - x_{i - 1}) (x - x_{i + 1}) \dots (x - x_{n})}{(x_{i} - x_{0}) \dots (x_{i} - x_{i - 1}) (x_{i} - x_{i + 1}) \dots (x_{i} - x_{n})} \\ (2) \\ (3) & = \underset{j \neq i}{\prod_{j = 0}^{n}} \frac{x - x_{j}}{x_{i} - x_{j}} (i = 0, 1, \dots, n) \end{aligned}

其中 $l_{i} (x)$ 是n次多项式，满足

l_{i} (x_{j}) = {\begin{cases} 0, j \neq i \\ 1, j = i \end{cases}

i,j相同时， $l_{i} (x)$ 式上下一致，结果为1

i,j不同时，$li(x) $中一定存在$ (x_j-x?)=0$，结果为0

拉格朗日插值函数：

L_{n} (x) = \sum_{i = 0}^{n} y_{i} l_{i} (x) = \sum_{i = 0}^{n} y_{i} (\underset{j \neq i}{\prod_{j = 0}^{n}} \frac{x - x_{j}}{x_{i} - x_{j}}) (i = 0, 1, \dots, n)

（三）样条插值

部分问题对所得插值函数的光滑性有较高的要求，不仅要连续，而且要有连续的曲率

样条（Spline）原指富有弹性的细木条或金属条。利用将已知点连接成一条光滑曲线（称为样条曲线），并使连接点处具有连续的曲率，三次样条插值就是如此抽象出来的

数学上将具有一定光滑性的分段多项式称为样条函数。具体来说，给定区间[a,b]的一个分划

Δ : a = x_{0} < x_{1} < \dots < x_{n - 1} < x_{n} = b

如果函数S(x)满足

在每一个小区间$[xi,x{i+1}]\ (i = 0,1,\cdots,n-1)$上S(x)是m次多项式
S(x)在 $[a, b]$ 上具有m-1阶连续导数

则称S(x)为关于分划 $Δ$ 的m次样条函数，其图形为m次样条曲线。显然，折现是依次样条曲线

三、MATLAB插值工具箱:notebook_with_decorative_cover:

（一）一维插值函数 interp1()

Code : <ValName> = interp1(x0, y0, x, <method>)

其中<method>指插值的方法，默认为线性插值(‘linear’)，其可选值有：

< method >	Description	Sp
‘linear’	线性插值	None
‘nearest’	最近项目插值	None
‘spline’	立方(三次)样条插值	保证函数平滑
‘cubic’	立方(三次)插值	保持函数凸凹性

其中所有的插值方法都要求x0是单调的

当x0为等距的时候可以使用快速插值法

method格式：'*linear' , '*nearest' , '*spline' , '*cubic'

（二）一维三次样条插值

CodeStyle01 : <ValName> = interp1(x0, y0, x, 'spline')

CodeStyle02 : <ValName> = spline(x0, y0, x)

Recommand : csape()

CodeStyle03 : pp = csape(x0, y0, <conds>)

<ValName> = fnval(pp, x)

CodeStyle04 : pp = csape(x0, y0, <conds>, <valconds>)

<ValName> = fnval(pp, x)

提倡使用函数csape，其返回值是pp形式，为了求出插值点的函数值，需要调用函数fnval

pp = csape(x0, y0)使用默认的边界条件（Lagrange边界条件）

<conds>为指定插值的边界条件，可选值如下

< conds > Description < valconds >

‘complete’ 边界为一阶导数需提供

‘not-a-knot’ 非扭结条件

‘periodic’ 周期条件

‘second’ 边界为二阶导数(缺省值为[0,0]) 需提供

‘variational’ 设置边界的二阶导数值为[0,0] 对于一些特殊的边界条件，可以通过conds的一个1 $\times$ 2矩阵来表示，conds元素的取值为0,1,2

Eg: conds = [2,1] 左边界是二阶导数，右边界是一阶导数，相应的值由valconds给出

<valconds>为<conds>方法所需的值，若不填写则按照缺省情况处理

< conds >	Description	< valconds >
‘complete’	边界为一阶导数	需提供
‘not-a-knot’	非扭结条件
‘periodic’	周期条件
‘second’	边界为二阶导数(缺省值为[0,0])	需提供
‘variational’	设置边界的二阶导数值为[0,0]	对于一些特殊的边界条件，可以通过conds的一个1 $\times$ 2矩阵来表示，conds元素的取值为0,1,2

Eg5.1 机床加工

（三）二维插值

1. 插值节点为网格节点

已知m $\times$ n个节点：$(xi,y_j,z{ij})\ (i = 1,2,\cdots,m;\ j=1,2,\cdots,n)$，

且 $x_{1} < \dots < x_{m}; y_{1} < \dots < y_{n}$

求点 $(x, y)$ 处的插值

（1）二维插值函数

Code : <ValName> = interp2(x0, y0, z0, x, y, <method>)

其中<method>指插值的方法，同一维的参数相同

x0,y0分别为m维和n维向量，表示节点，z0为n $\times$ m维度矩阵，表示节点的值，x,y为一维数组，指插值点，

x与y应当是方向不同的向量，即一个是行向量，一个是列向量（y需要转置）

输出内容为矩阵，Row=y的维数，Col=x的维数，表示得到的插值

（2）二维三次样条插值

Code : pp = csape({x0, y0}, z0, conds, valconds)

<ValName> = fnval(pp, {x, y})

对其内容的解释同上

（3）显示相应的三维网格图

Code : mesh(y,x,z)

mesh函数的x为列，y为行，需要将xy反写

Eg5.3 丘陵曲面最高点

2. 插值节点为散乱节点

已知n个节点 $(x_{i}, y_{i}, z_{i}) (i = 1, 2, \dots, n)$ ，求点 $(x, y)$ 处的插值z

Code：<ValName> = griddata(x, y, z, XI, YI)

其中x、y、z均为n维向量，指明所给出数据的横坐标、纵坐标和竖坐标

向量XI、YI是给定的网格点的横坐标和纵坐标

返回值<ValName>为网格(XI, YI)处的函数值，但是行为YI，列为XI，进行mesh的时候不需要反转XY

XI与YI应当是方向不同的向量，一行一列（YI需要转置）

MATLAB中插值时，外插值是不确定的，可以使用混合插值的方法，把不确定的插值换成最近点的插值结果
matlab
1
2
3
4
ZI1=griddata(x,y,z,XI,YI','cubic'); %立方插值
ZI2=griddata(x,y,z,XI,YI','nearest'); %最近点插值
ZI = ZI1; %
ZI(isnan(ZI1)) = ZI2(isnan(ZI1));

四、实战例题

Eg5.1 机床加工

待加工零件的外形根据工艺要求由一组数据 $(x, y)$ 给出（在平面情况下），用程控铣床加工时每一刀只能沿 $x$ 方向和 $y$ 方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的 $(x, y)$ 坐标。
Data中给出的 $x, y$ 数据位于机翼断面的下轮廓线上，假设需要得到 $x$ 坐标==每改变0.1时==的 $y$ 坐标。试完成加工所需数据，==画出曲线==，并==求出 $x = 0$ 处的曲线斜率和 $13 \leq x \leq 15$ 范围内y的最小值==。要求用==分段线性==和==三次样条==两种插值方法计算。

Data：

x	0	3	5	7	9	11	12	13	14	15
y	0	1.2	1.7	2.0	2.1	2.0	1.8	1.2	1.0	1.6

matlab

clear
clc
%% Data:
x0 = [0 3 5 7 9 11 12 13 14 15];
y0 = [0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
x = [min(x0):0.1:max(x0)];
%% Linear
y_linear = interp1(x0, y0, x);
%% Spline
y_spline = interp1(x0, y0, x, 'spline');
pp_lagrange = csape(x0, y0);
y_spline_lagrange = fnval(pp_lagrange, x);
pp_second = csape(x0, y0, 'second');
y_spline_second = fnval(pp_second, x);
%% DispFigure
% linear
figure('Position', [800,100,1000,800]);
subplot(2,2,1)
plot(x,y_linear, 'r',x0, y0,'bp')
title('Interp1-Linear')
% spline
subplot(2,2,2)
plot(x,y_spline, 'r',x0, y0,'bp')
title('Interp1-Spline')
subplot(2,2,3)
plot(x,y_spline_lagrange, 'r',x0, y0,'bp')
title('Csape-Lagrange')
subplot(2,2,4)
plot(x,y_spline_second, 'r',x0, y0,'bp')
title('Csape-Second')
%% Other Tasks
% 求x=0处的斜率
dx = diff(x); % diff逐差函数
dy = diff(y_spline_lagrange);
dy_dx = dy./dx;
dy_dx_0 = dy_dx(1);
% 求13~15范围内y的最小值
ymin = min(y_spline_lagrange(131:151)); % 找出ymin
yloc = find(y_spline_lagrange == ymin); % find ymin的index
xloc = x(yloc); % 用ymin的index定位xloc
disp(['k(0) = ' num2str(dy_dx_0)])
disp(['xloc: ' num2str(xloc)])
disp(['ymin: ' num2str(ymin)])

output

plaintext

1
2
3

k(0) = 0.49716
xloc: 13.8
ymin: 0.98511

Pic

Eg5.2 求速度曲线的位移

已知速度曲线 $v (t)$ 上的四个数据点如下表，用三次样条插值求位移
$S = \int_{0.15}^{0.18} v (t) d t$

Data:

$t$	0.15	0.16	0.17	0.18
$v (t)$	3.5	1.5	2.5	2.8

matlab

clear
clc
%% DataLoad
t0 = 0.15:0.01:0.18;
v0 = [3.5,1.5,2.5,2.8];
%% Spline
pp = csape(t0, v0); % 使用默认边界条件Lagrange
format long g % 设置长小数格式
xs = pp.coefs % 显示每个区间上三次多项式的系数
s = integral(@(t)ppval(pp,t),0.15,0.18) % 求数值积分
    % @(t) 为匿名函数的书写格式 
    % integral(fx,xmin,ymin) 求数值积分
        % 使用全局自适应积分和默认误差容限在 xmin 至 xmax 间以数值形式为函数
    % ppval(pp,xq)计算分段多项式:查询点 xq 处计算分段多项式 pp
format % 恢复短小数小时格式

output

plaintext

xs =     -616666.666666667                     33500         -473.333333333334                       3.5
         -616666.666666667                     15000          11.6666666666671                       1.5
         -616666.666666668         -3499.99999999999          126.666666666667                       2.5
s =0.068625

Eg5.3 丘陵曲面最高点

在一丘陵地带测量高程，==x和y方向每隔100米==测一个点，得高程如表，试插值一==曲面==，确定合适的模型，并由此找出==最高点和该点的高程==

Data:

y\x	100	200	300	400	500
100	636	697	624	478	450
200	698	712	630	478	420
300	680	674	598	412	400
400	662	626	552	334	310

注意此处直接列写矩阵，则xy与正常的对应相反

matlab

clear
clc
%% DataLoad
x0 = 100:100:500;
y0 = 100:100:400;
z0 = [636 697 624 478 450
      698 712 630 478 420
      680 674 598 412 400
      662 626 552 334 310]';
% 此处的z0两坐标方向交叉，需要转置
step = 5;
x = 100:step:500;
y = 100:step:400;
%% interp
% interp2要求矩阵和向量交叉对应，且y为列向量
z1 = interp2(x0,y0,z0',x,y');
z1 = z1'; % 转置z1回正常的对应
% csape正常进行
pp = csape({x0,y0},z0);
z2 = fnval(pp,{x,y});
%% DispMesh
% mesh函数的x为列，y为行，需要将xy反写
figure('Position',[500,300,1000,600])
subplot(1,2,1)
mesh(y,x,z1) 
title('Interp1-Linear','FontSize',20)
subplot(1,2,2)
mesh(y,x,z2)
title('Spline','FontSize',20)
%% MAXValue
% interp-Linear
[maxls1,xlocls1] = max(z1); % 按列取max
[vmax1,yloc1] = max(maxls1); % 按行取max
info1 = ["interp-Linear",x(xlocls1(yloc1)),y(yloc1),vmax1]; 
    % 行loc索引列ls的loc,并用loc索引相应的x,y的值
% Spline
[maxls2,xlocls2] = max(z2);
[vmax2,yloc2] = max(maxls2);
info2 = ["Spline",x(xlocls2(yloc2)),y(yloc2),vmax2]; 
info = [info1;info2];
% 更好的位置索引方法
[i,j]=find(z2==max(max(z2)));  %找最高点的地址
xm=x(i),ym=y(j),zm=z2(i,j);  %求最高点的坐标 
[zm,xm,ym]
%% DispMAX
disp(['Step:' num2str(step)])
for i = 1:2
    disp(['Type:' num2str(info(i,1))]) 
    disp(['max = ' num2str(info(i,4)) ...
        '    x = ' num2str(info(i,2)) ...
        '    y = ' num2str(info(i,3))])
end

output

plaintext

  720.7217  165.0000  175.0000
Step:5
Type:interp-Linear
max = 712    x = 200    y = 200
Type:Spline
max = 720.7217    x = 165    y = 175

Pic

Eg5.4 海底图像

在某海域测得一些点（x,y）处的水深z由表5.4给出，在适当的矩形区域内画出海底曲面的图形。

Data：

x=[129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5];

y=[7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5];

z=-[4,8,6,8,6,8,8,9,9,8,8,9,4,9];

matlab

clear
clc
%% DataLoad
x=[129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5];
y=[7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5];
z=-[4,8,6,8,6,8,8,9,9,8,8,9,4,9];
step = 1;
XI = min(x):step:max(x)
YI = min(y):step:max(y)
%% grid
ZI1=griddata(x,y,z,XI,YI','cubic'); %立方插值
ZI2=griddata(x,y,z,XI,YI','nearest'); %最近点插值
ZI = ZI1; % 设置ZI初始值为ZI1
ZI(isnan(ZI1)) = ZI2(isnan(ZI1)); % 用ZI2填充ZI1外插值(为NaN)的部分
figure
mesh(XI,YI,ZI)
hold on
plot3(x,y,z,'or')

Pic

TestDoc

插值

一、插值定义

二、常见插值方法

（一）分段线性插值

（二）拉格朗日插值多项式

（三）样条插值

三、MATLAB插值工具箱:notebook_with_decorative_cover:

（一）一维插值函数 interp1()

（二）一维三次样条插值

（三）二维插值

1. 插值节点为网格节点

（1）二维插值函数

（2）二维三次样条插值

（3）显示相应的三维网格图

2. 插值节点为散乱节点

四、实战例题

Eg5.1 机床加工

Eg5.2 求速度曲线的位移

Eg5.3 丘陵曲面最高点

Eg5.4 海底图像

微信扫一扫：分享