所属学科
计算方法、数值分析、数值代数、数值逼近
求一个方程、积分、微分方程、方程组等的数值解
问题
美国的人口普查每10年举行一次。下表列出来了从1940年到1990年的人口（按千人计）。
请估计1965年的人口。
工程实践和科学实验中，常常需要从一组实验观察数据(i=1,2,…，n)中，求自变量x与因变量y的一个近似的函数关系式y=f(x)。
例如观测行星的运动，只能得到某时刻ti所对应的行星位置si（用经纬度表示），想知道行星在任何时刻t的位置。
这类问题都需要我们通过部分已有的数据去估计未知的数据，去实现数据的数值逼近，其中两个方法就是插值和拟合。
区别
首先插值和拟合都是根据某个未知函数（或已知但难于求解的函数）的几个已知数据点求出变化规律和特征相似的近似曲线的过程。但是插值法要求的是近似的曲线需要完全经过数据点，而拟合则是得到最接近的结果，强调最小方差的概念。
插值：满足 即已有的数据都在函数上
拟合：满足 ⇒ 即该曲线在某种准则下与所有的数据点最为接近，即曲线拟合的最好（最小化损失函数）
插值
代数多项式简单便于计算，所以用代数多项式近视地表示满足n个点的函数关系式y=f(x)。
拉格朗日插值
n次拉格朗日插值公式：
%拉格朗日插值
function yh=lagrange (x,y,xh)
n = length(x);
m = length(xh);
yh = zeros(1,m);
c1 = ones(n-1,1);
c2 = ones(1,m);
for i=1:n
xp = x([1:i-1 i+1:n]);
yh = yh + y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));
end牛顿均差插值
Runge现象和分段低次插值
clear;clc;
hold off
x=[-5:1:5];
y=1./(1+x.^2);
x0=[-5:0.1:5];
y0=lagrange(x,y,x0);
y1=1./(1+x0.^2);
plot(x0,y0)
hold on
plot(x0,y1,'b:')
legend('插值曲线','原数据曲线')
MATLAB实现
一维数据插值（一元函数）
yi=interp1(x,y,xi,method)
%1.'linear':分段线性插值
%2.'pchip':分段三次多项式插值
%3.'spline':三次样条插值
%4.'nearest':最近邻近区域插值（阶梯函数）
例题
hold off
xx=1:1:5;
yx=[3.5,4.6,5.5,3.2,2];
xxi=1:0.5:5;
fo=interp1(xx,yx,xxi);
f1=interp1(xx,yx,xxi,'linear');
f2=interp1(xx,yx,xxi,'pchip');
f3=interp1(xx,yx,xxi,'spline');
f4=interp1(xx,yx,xxi,'nearest');
f0,f1,f2,f3,f4三次样条插值
yi=spline(x,y,xi) 求某个点的数据
例题
xx=1:1:5;
yx=[3.5,4.6,5.5,3.2,2];
xxi=1:0.5:5;
f0=interp1(xx,yx,xxi)
f1=interp1(xx,yx,xxi,'spline')
yi=spline(xx,yx,xxi)
pp=spline(x,y) 插值函数的数据
例题
xx=1:1:5;
yx=[3.5,4.6,5.5,3.2,2];
xxi=1:0.5:5;
f0=interp1(xx,yx,xxi)
ps=spline(xx,yx);
yyi=ppval(ps,xxi)二维数据插值（二元函数）
曲面插值
zi=interp2(x,y,z,xi,yi,method)
例题
%某实验对一根长10m的钢轨进行热源的温度传播测试。用x表示测试点0：2.5：10（m），
%用h表示测试时间0：30：60（s），用T表示测试所得个点的温度
%试用双线性插值在一分钟内每隔20s、钢轨每m处的温度Ti
x=0:2.5:10;
h=[0:30:60]';
T=[95 14 0 0 0;88 48 32 12 6;67 64 54 48 41];
mesh(x,h,T);
xi=[0:10];
hi=[0:2:60]';
Ti=interp2(x,h,T,xi,hi);
mesh(xi,hi,Ti)
插值基点为散乱节点
cz=griddata(x,y,z,cx,cy,method)
%cx,cy一个为行向量一个为列向量
例题
%海底曲面图形
%在某海域测得一些点（x，y）处的水深z，在矩形区域（75，200）*（-50，150）内画出海底曲面图形）
clear;clc;
x=[129 140 103.5 88 185.5 195 105 158 108 77 81 162 162 118];
y=[7.5 142 23 147 23 138 86 -6.5 -81 3 57 -66 84 -34];
z=[-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9];
cx=linspace(75,200,100);
cy=linspace(-50,150,100);
cz=griddata(x,y,z,cx,cy','cubic')
meshz(cx,cy,cz)
contour3(cx,cy,cz,16)题集（可以将程序运行看看插值的效果）
1、
clear;clc;
x=[0.4,0.5,0.6,0.7,0.8];
y=[-0.916291,-0.693147,-0.510826,-0.356675,-0.223144];
X=[0.55,0.66];
Y=interp1(x,y,X)
Y1=log(X)
2、
clear;clc;
x = linspace(0, 2*pi, 12);
y = sin(x)
x1 = linspace(0, 2*pi, 30);
y1 = interp1(x,y,x1,'spline');
x2 = x1;
y2 = sin(x2);
plot(x,y,'r*')
hold on
plot(x1,y1,'b-')
plot(x2,y2,'k--')
legend('插值点','样条插值','解析解')
3、
clear;clc;
x=0:5;
y=1-cos(x).*exp(-x);
plot(x,y,'*')
X1=linspace(0,5,30);
X2=X1;
Y1=interp1(x,y,X1,'linear');
Y2=interp1(x,y,X2,'cubic');
X=X1;
Y=1-cos(X).*exp(-X);
hold on
plot(X1,Y1,'b-')
plot(X2,Y2,'m-.')
plot(X,Y,'k')
legend('插值点','线性插值','三次插值','解析解')
4、
clear;clc;
t = linspace(0, 2*pi, 36);
x = cos(t);
y = sin(t);
z = abs(cos(2*t));
t1 = linspace(0, 2*pi, 360);
x1 = interp1(t,x,t1);
y1 = interp1(t,y,t1);
z1 = interp1(t,z,t1);
x2 = cos(t1);
y2 = sin(t1);
z2 = abs(cos(2*t1));
plot3(x,y,z,'r*')
hold on
plot3(x1,y1,z1,'b-')
plot3(x2,y2,z2,'k--')
legend('插值点','线性插值','解析解')
box on
view(-63,38)
5、
clear;clc;
x=-5:5;
y=-5:5;
[X,Y]=meshgrid(x,y); %原始插值点
Z=exp(-((X-1).^2+(Y-2).^2))-exp(-((X+2).^2+(Y+1).^2)); %原始插值点
surf(X,Y,Z) %绘制原始插值点曲面
title('插值点曲面')
x1=linspace(-5,5,50);
y1=linspace(-5,5,50);
[X1,Y1]=meshgrid(x1,y1);
Z1=interp2(X,Y,Z,X1,Y1,'pchip'); %三次插值运算
figure(2)
surf(X1,Y1,Z1)
title('三次插值曲面') %绘制三次插值曲面
figure(3)
[XX,YY]=meshgrid(x1,y1);
ZZ=exp(-((XX-1).^2+(YY-2).^2))-exp(-((XX+2).^2+(YY+1).^2));
surf(XX,YY,ZZ)
title('解析解曲面') %绘制解析解插值曲面
6、
clear;clc;
x = linspace(0, 2*pi, 6);
y = linspace(0, 2*pi, 6);
[X,Y] = meshgrid(x,y);
Z = cos(X)+sin(Y);
surf(X,Y,Z)
x1 = linspace(0, 2*pi, 36);
y1 = linspace(0, 2*pi, 36);
[X1,Y1] = meshgrid(x1,y1);
Z1 = interp2(X,Y,Z,X1,Y1,'cubic');
figure (2)
surf(X1,Y1,Z1)
7、
clear;clc;
x=-5:1:5;
y=exp(-abs(x));
N=length(x);
X=linspace(-5,5,11*N);
Y=spline(x,y,X);
plot(x,y,'*')
hold on
plot(X,Y,'r--')
X1=X;
Y1=exp(-abs(X1));
plot(X1,Y1)
legend('插值点','样条插值','解析解')
8、
clear; clc;
t = linspace(0, 2*pi, 9);
x = cos(t);
y = sin(t);
z = t.^2;
t1 = linspace(0, 2*pi, 36);
x1 = interp1(t, x, t1);
y1 = interp1(t, y, t1);
z1 = interp1(t, z, t1);
t2 = t1;
x2 = spline(t, x, t2);
y2 = spline(t, y, t2);
z2 = spline(t, z, t2);
plot3(x, y, z, 'r*');
hold on
plot3(x1,y1,z1, 'm-.')
plot3(x2, y2, z2, 'b-')
box on
legend('插值点','线性插值','样条插值')
view(-50,26)拟合
已知一组（二维）数据，即平面上的n个点互不相同。寻求一个函数（曲线），使在某种准则下与所有数据点最为接近，即曲线拟合得最好。其中线性最小二乘法是解决曲线拟合最常用的方法。
MATLAB实现
polyfit(x,y,N) %采用最小二乘法构造一个N次多项式
P=polyfit(x,y,n) %用n阶多项式拟合x，y向量给定的数据
PA=polyval(p,xi) %求xi点上的拟合函数的近视值
px=poly2str(p,'x')
例题
x=[0.5 1.0 1.5 2.0 2.5 3.0];
y=[1.75 2.45 3.81 4.80 7.00 8.60];
a1=polyfit(x,y,1)
a2=polyfit(x,y,2)
a3=polyfit(x,y,3)
x1=[0.5:0.05:3.0];
y1=a1(2)+a1(1)*x1;
y2=a2(3)+a2(2)*x1+a2(1)*x1.^2;
y3=a3(4)+a3(3)*x1+a3(2)*x1.^2+a3(1)*x1.^3;
plot(x,y,'*')
hold on
plot(x1,y1,'b:',x1,y2,'k',x1,y3,'g')
legend('原数据图','一次拟合','二次拟合','三次拟合')
%非线性最小二乘曲线拟合 lsqnonline lsqcurvefit
[x,resnorm]=lsqcurvefit(fun,x0,xdata,ydata)
%fun:需要拟合的函数
%x0:对函数中参数的初始值
%xdata:已知x轴的数据
%ydata:已知y轴的数据
%输出
%x:最优解
%resnorm：误差的平方和
lsqnonline()
例题
% y=ax+bx^2*e^(-cx)+d
clear;clc;
f=inline('a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4)','a','x')
x=0.1:0.1:1;
y=[2.3201 2.6470 2.9707 3.2885 3.6008 3.9090 4.2147 4.5191 4.8232 5.1275];
[xx,res]=lsqcurvefit(f,ones(1,4),x,y)题集
1、
clear;clc;
x=[0.5, 1.0, 1.5, 2.0, 2.5, 3.0];
y=[1.75, 2.45, 3.81, 4.80, 7.00, 8.60];
a=polyfit(x, y, 2);
x1=linspace(0.5, 3.0, 20);
y1=a(1)*x1.^2+a(2)*x1+a(3);
plot(x,y,'*')
hold on
plot(x1,y1,'r-')
legend('数据点','拟合曲线')
2、
clear;clc;
x = linspace(0, pi/2, 6);
y = sin(x);
a = polyfit(x,y,2);
x1 = linspace(0, pi/2, 72);
y1 = a(1)*x1.^2+a(2)*x1+a(3);
b = polyfit(x,y,3);
x2 = x1;
y2 = b(1)*x1.^3+b(2)*x1.^2+b(3)*x1+b(4);
plot(x,y,'r*')
hold on
plot(x1,y1,'b-')
plot(x2,y2,'k--')
legend('插值点','二次多项式','三次多项式')作业：matlab中插值和拟合函数简单的使用
1. 已知1950年到1990年间每隔10年，服务年限从10年到30年每隔10年的劳动报酬表如下：
表：某企业工作人员的月平均工资（元）
服务年限
年份 10 20 30
1950 150.697 169.592 187.652
1960 179.323 195.072 250.287
1970 203.212 239.092 322.767
1980 226.505 273.706 426.730
1990 249.633 370.281 598.243
试计算1975年时，15年工龄的工作人员平均工资。
2. 根据以下数表，求拟合函数y=ax+b.
x 1 2 3 4 5 6 7
y 5.126 8.119 11.4984 14.9597 17.3404 20.5853 23.2238