正文
Matlab随笔之插值与拟合（下）
小程序：扫一扫查出行【扫一扫了解最新限行尾号】
复制小程序
原文:Matlab随笔之插值与拟合（下）
1、二维插值之插值节点为网格节点
已知m x n个节点：（xi，yj，zij）(i=1…m,j=1…n)，且xi,yi递增。求（x，y）处的插值z。
Matlab可以直接调用interp2(x0,y0,z0,x,y,`method`)
其中 x0,y0 分别为 m 维和 n 维向量，表示节点， z0 为 n × m 维矩阵，表示节点值， x,y
为一维数组，表示插值点， x 与 y 应是方向不同的向量，即一个是行向量，另一个是列
向量， z 为矩阵，它的行数为 y 的维数，列数为 x 的维数，表示得到的插值， 'method'
的用法同上面的一维插值。
如果是三次样条插值，可以使用命令
pp=csape({x0,y0},z0,conds,valconds), z=fnval(pp,{x,y})
其中 x0,y0 分别为 m 维和 n 维向量， z0 为 m × n 维矩阵， z 为矩阵，它的行数为 x 的维
数，列数为 y 的维数，表示得到的插值，具体使用方法同一维插值。
eg:
(1)、用interp2函数插值：
x=::;
y=::;
z=[ ];
p=::;
q=::;
q=q';%须为列向量
z0=interp2(x,y,z,p,q);%分段线性插值
z1=interp2(x,y,z,p,q,'spline');%三次线条插值
subplot(,,);
mesh(p,q,z0);
title('分段线性插值');
subplot(,,);
mesh(p,q,z1);
title('三次线条插值');
%可以观察出，三次线条插值的图像更平滑
运行结果：
(2)、用csape函数插值：
x=::;
y=::;
z=[ ];
p=::;
q=::;
q=q';
%三次线条插值
pp=csape({x,y},z');%注意跟interp2的区别,有个转置
z0=fnval(pp,{p,q});
mesh(p,q,z0');%注意跟interp2的区别,有个转置
title('三次线条插值');
运行结果：
2、二维插值之插值节点为散乱节点
已知 n 个节点： ( xi , yi , zi )(i = 1,2,…, n) ，求点 (x, y) 处的插值 z 。
对上述问题， Matlab 中提供了插值函数 griddata，其格式为：
ZI = GRIDDATA(X,Y,Z,XI,YI)
其中 X、 Y、 Z 均为 n 维向量，指明所给数据点的横坐标、纵坐标和竖坐标。向量 XI、
YI 是给定的网格点的横坐标和纵坐标，返回值 ZI 为网格（ XI， YI）处的函数值。 XI
与 YI 应是方向不同的向量，即一个是行向量，另一个是列向量。
eg:
%散乱节点的二维插值
x=[ 103.5 185.5 157.5 107.5 117.5];
y=[7.5 141.5 22.5 137.5 85.5 -6.5 - 56.5 -66.5 -33.5];
z=-[ ];
x0=[::];
y0=[-::]';
z0=griddata(x,y,z,x0,y0,'cubic');%保凹凸性3次插值
%[xx,yy]=meshgrid(x0,y0);无需采样，故不需要该函数
mesh(x0,y0,z0);
运行结果:
在上述问题中,补上寻找最大值的程序:
%max(z0)返回一个行向量，向量的第i个元素是矩阵A的第i列上的最大值
%find(A) 寻找矩阵A非零元素下标,返回矩阵A中非零元素所在位置
%[i,j,v]=find(A)返回矩阵A中非零元素所在的行i,列j,和元素的值v(按所在位置先后顺序输出)
[p,q]=find(z0==max(max(z0)));
zmax=z0(p,q)
3、最小二乘法实现曲线拟合（1）用最小二乘法求一个形如 y = a + bx^ 2 的经验公式：
%等价于[,x^][a;b]=y,转换成解超定方程问题，超定方程的解是根据最小二乘法得来的x=[ ]';
y=[19.0 32.3 49.0 73.3 97.8]';
r=[ones(,),x.^]
ab=r\y
x0=:0.1:;
y0=ab()+ab()*x0.^;
plot(x,y,'o',x0,y0,'r')
（2）多项式拟合
%等价于[,x^][a;b]=y,转换成解超定方程问题，超定方程的解是根据最小二乘法得来的x=[ ]';
y=[19.0 32.3 49.0 73.3 97.8]';
r=[ones(,),x.^]
ab=r\y
x0=:0.1:;
y0=ab()+ab()*x0.^;
plot(x,y,'o',x0,y0,'r')
%a=polyfit(x,y,n)用多项式求过已知点的表达式，其中x为源数据点对应的横坐标，可为行向量、矩阵，y为源数据点对应的纵坐标，可为行向量、矩阵，n为你要拟合的阶数，一阶直线拟合，二阶抛物线拟合，并非阶次越高越好，看拟合情况而定，a为m+1的行向量。polyfit函数的数学基础是最小二乘法曲线拟合原理，所得到的函数值在基点处的值与原来点的坐标偏差最小，常用于数据拟合，polyfit 做出来的值从左到右表示从高次到低次的多项式系数。
如果要求拟合函数在x`点的函数值，可以调用polyval(a,x`)函数
eg：
x0=[ ];
y0=[ ];
%画出散点图
plot(x0,y0,'ro');
hold on
%用线性拟合
p=polyfit(x0,y0,);
z0=polyval(p,x0);
plot(x0,z0);
运行结果：
4、最小二乘优化 (最小二乘:least square)
1、lsqlin函数
eg：
%拟合形如y=a+bx^2的函数
%采样点
x=[ ]';
y=[ 32.3 73.3 97.8]';
r=[ones(,),x.^];
ab=lsqlin(r,y)
x0=:0.1:;
y0=ab()+ab()*x0.^;
plot(x,y,'o',x0,y0,'r')
运行结果：
5、曲线拟合与函数逼近
eg：
求 f(x) =cos x, （-pi/2<=x<=pi/2） 在 H = Span{1, x^2 , x^4} 中的最佳平方逼近多项式。
程序如下：
syms x%定义符号数值
base=[,x^,x^];
y1=base.'*base
y2=cos(x)*base.'
r1=int(y1,-pi/,pi/)
r2=int(y2,-pi/,pi/)
a=r1\r2%a为符号数值
xishu1=double(a)%化简符号数值
digits()%设置符号数值的精度
xishu2=vpa(a)%任意精度（符号类）数值
运行结果：
y1 = [ 1, x^2, x^4]
[ x^2, x^4, x^6]
[ x^4, x^6, x^8] y2 = cos(x)
x^2*cos(x)
x^4*cos(x) r1 = [ pi, pi^3/12, pi^5/80]
[ pi^3/12, pi^5/80, pi^7/448]
[ pi^5/80, pi^7/448, pi^9/2304] r2 = 2
pi^2/2 - 4
pi^4/8 - 6*pi^2 + 48 a = (15*(pi^4 - 308*pi^2 + 3024))/(4*pi^5)
-(210*(pi^4 - 228*pi^2 + 2160))/pi^7
(1260*(pi^4 - 180*pi^2 + 1680))/pi^9 xishu1 = 0.9996
-0.4964
0.0372xishu2 = 0.99957952
-0.49639233
0.037209327 >>
所以y的最佳平方逼近多项式为y=0.9996-0.4964x^2+0.0372x^4
