TestDoc

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

插值

[toc]

一、插值定义

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

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

二、常见插值方法

已知$n+1$个点$(x_i,y_i)\ (i=0,1,…,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)$，其中

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

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

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

拉格朗日插值的基函数为

其中$l_i(x)$是n次多项式，满足

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

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

拉格朗日插值函数：

（三）样条插值

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

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

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

如果函数S(x)满足

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

则称S(x)为关于分划$\Delta$的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>方法所需的值，若不填写则按照缺省情况处理

Eg5.1 机床加工

（三）二维插值

1. 插值节点为网格节点

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

且$x_1<\cdots<x_m;\ y_1<\cdots<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,\cdots,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中插值时，外插值是不确定的，可以使用混合插值的方法，把不确定的插值换成最近点的插值结果

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\le x\le15$范围内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 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
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

1
2
3
k(0) = 0.49716
xloc: 13.8
ymin: 0.98511

Pic

Eg5.2 求速度曲线的位移

已知速度曲线$v(t)$上的四个数据点如下表，用三次样条插值求位移

Data:

$t$ 0.15 0.16 0.17 0.18
$v(t)$ 3.5 1.5 2.5 2.8 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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

1
2
3
4
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与正常的对应相反

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
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

1
2
3
4
5
6
  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];

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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