matlab中函数拟合方法—个人总结.doc

目 录
一维插值方案 2
二维数据内插值(表格查找) 3
等高线 4
三维曲面 5
等高线2 6
三维曲面2 7
matlab绘制温度场(尚未深入研究) 13
二维曲线（非线性）拟合步骤 18
三维曲线（非线性）拟合步骤 19
三维曲线的画法 20
三维曲面的画法 21
画三维图3 只有点的数据,没有函数关系式 23
空间点拟合的基本原理 27
空间点拟合的最小二乘法 28
曲面生成后再进行多项式拟合 37
六点生成曲面 38
四点生成平面 39
用三维离散点拟合光滑曲面1 40
用三维离散点拟合光滑曲面2 40

一维插值方案
clear
year = 1900:10:2010;
product = [75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 256.344 267.893 ]
p1995 = interp1(year,product,1995)
%使用一维数据内插值(该题中只能在1900和2010之间进行插值,大于2010和小于1900都%无效)命令
x = 1900:1:2010
y = interp1(year,product,x,'spine');
plot(year,product,'o',x,y)

插值说明:
interp1(x,Y,xi,method) %用指定的算法计算插值：
’nearest’：最近邻点插值，直接完成计算；
’linear’：线性插值（缺省方式），直接完成计算；
’spine’：三次样条函数插值。对于该方法，命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函
数。命令spline 用它们执行三次样条函数插值；
’pchip’：分段三次Hermite 插值。对于该方法，命令interp1 调用函数pchip，用于对向量x 与y 执行分段三次内插值。该方法保留单调性与
数据的外形；
’cubic’：与’pchip’操作相同；
’v5cubic’：在MATLAB 5.0 中的三次插值。
对于超出x 范围的xi 的分量，使用方法’nearest’、’linear’、’v5cubic’的插值算法，相应地将返回NaN。对其他的方法，interp1 将对超出的分量执行外插值算法。
yi = interp1(x,Y,xi,method,'extrap') %对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。
yi = interp1(x,Y,xi,method,extrapval) %确定超出x 范围的xi 中的分量的外插值extrapval，其值通常取NaN 或0。
例1

clear;
x = 0:10; y = x.*sin(x);
xx = 0:.25:10; yy = interp1(x,y,xx)
plot(x,y,'kd',xx,yy)

interp2
二维数据内插值(表格查找)
[X,Y] = meshgrid(-3:.25:3);
Z = peaks(X,Y);
[XI,YI] = meshgrid(-3:.125:3);
ZZ = interp2(X,Y,Z,XI,YI);
surfl(X,Y,Z);hold on;
surfl(XI,YI,ZZ+15)
axis([-3 3 -3 3 -5 20]);
shading flat
hold off

功能 三维数据插值interp3（查表）
[x,y,z,v] = flow(20);
[xx,yy,zz] = meshgrid(.1:.25:10, -3:.25:3, -3:.25:3);
vv = interp3(x,y,z,v,xx,yy,zz);
slice(xx,yy,zz,vv,[6 9.5],[1 2],[-2 .2]); shading interp;colormap cool

等高线
clear
Z=peaks
for w=1:1:100
V=[w/10,0,w/10]
contour(Z,V)
%C=contour(Z,V)
%Clabel(C)
Hold on
title('等高线及其标注')
end
end

三维曲面
x=0:10
y=0:.1:1
[d,B]=meshgrid(x,y)
z=1./(B.*d.^2+1);
surf(B,d,z)

x=0:0.05:10
y=0:0.05:1
[X,Y]=meshgrid(x,y)
Z=( X.^3+ 3.*Y.^2+5*Y); %Z=( X.^2+ 3.*Y.^3+5*Y);%
surf(X,Y,Z)
%一张普通的三维曲面,有时需要旋转一下才能看到下图的结果;

x=0:0.05:1
y=0:0.05:1
[X,Y]=meshgrid(x,y)
Z=( X.^2-Y.^2);% Z=( 4*X.^3*Y-4*X.*Y.^3);
surf(X,Y,Z) %一张普通的三维曲面,有时需要旋转一下才能看到下图的结果;

等高线2
clear
x=-2:0.1:2
y=-2:0.1:2
[X,Y]=meshgrid(x,y)
Z=(X.^2+Y.^2).^0.5
for w=1:1:100
V=[w/3,w/pi,w/3]
contour(Z,V)
hold on
end

三维曲面2

clear
x=-5:0.05:5
y=-5:0.05:5
[X,Y]=meshgrid(x,y)
Z=1./((X+1).^2+(Y+1).^2+1)-1.5./((X-1).^2+(Y-1).^2+1)
mesh(X,Y,Z)

clear;
A=[1.486,3.059,0.1;2.121,4.041,0.1;2.570,3.959,0.1;3.439,4.396,0.1;
4.505,3.012,0.1;3.402,1.604,0.1;2.570,2.065,0.1;2.150,1.970,0.1;
1.794,3.059,0.2;2.121,3.615,0.2;2.570,3.473,0.2;3.421,4.160,0.2;
4.271,3.036,0.2;3.411,1.876,0.2;2.561,2.562,0.2;2.179,2.420,0.2;
2.757,3.024,0.3;3.439,3.970,0.3;4.084,3.036,0.3;3.402,2.077,0.3;
2.879,3.036,0.4;3.421,3.793,0.4;3.953,3.036,0.4;3.402,2.219,0.4;
3.000,3.047,0.5;3.430,3.639,0.5;3.822,3.012,0.5;3.411,2.385,0.5;
3.103,3.012,0.6;3.430,3.462,0.6;3.710,3.036,0.6;3.402,2.562,0.6;
3.224,3.047,0.7;3.411,3.260,0.7;3.542,3.024,0.7;3.393,2.763,0.7];
x=A(:,1);y=A(:,2);z=A(:,3);
scatter(x,y,5,z)%散点图
figure
[X,Y,Z]=griddata(x,y,z,linspace(1.486,4.271)',linspace(1.604,4.276),'v4');%插值
pcolor(X,Y,Z);shading interp%伪彩色图
figure,contourf(X,Y,Z) %等高线图

clear;
A=[1.486,3.059,1858;2.121,4.041, 1858;2.570,3.959, 1858;3.439,4.396, 1858;
4.505,3.012, 1858;3.402,1.604, 1858;2.570,2.065, 1858;2.150,1.970, 1858;
1.794,3.059,2350;2.121,3.615, 2350;2.570,3.473, 2350;3.421,4.160, 2350;
4.271,3.036, 2350;3.411,1.876, 2350;2.561,2.562, 2350;2.179,2.420, 2350;
2.757,3.024, 2600;3.439,3.970, 2600;4.084,3.036, 2600;3.402,2.077, 2600;
2.879,3.036, 2849;3.421,3.793, 2849;3.953,3.036, 2849;3.402,2.219, 2849;
3.000,3.047, 3010;3.430,3.639, 3010;3.822,3.012, 3010;3.411,2.385, 3010;
3.103,3.012, 3345;3.430,3.462, 3345;3.710,3.036, 3345;3.402,2.562, 3345;
3.224,3.047, 3629;3.411,3.260, 3629;3.542,3.024, 3629;3.393,2.763, 3629];
x=A(:,1);y=A(:,2);z=A(:,3);
scatter(x,y,5,z)%散点图,5是点的大小
figure %打开显示图的界面
[X,Y,Z]=griddata(x,y,z,linspace(1.486,4.271)',linspace(1.604,4.276),'v4');%插值
pcolor(X,Y,Z);shading interp%伪彩色图
figure;contourf(X,Y,Z) %等高线图
figure;mesh(X,Y,Z)

A=[1.109,1.059,1718;2.021,0.841, 1758;2.870,0.359, 1858;4.039,0.196, 1838;
4.505,3.012, 3345;3.402,1.604, 3347;2.570,2.065, 3629;2.150,1.970, 3330;
1.794,3.059,2250;2.121,3.615, 3027;2.570,3.473, 2935;3.421,4.160, 1930;
4.271,3.036, 2050;3.411,1.876, 3144;2.561,2.562, 3739;2.179,2.420, 1950;
2.757,3.024, 3530;3.439,3.970, 2720;4.084,3.036, 2610;3.402,2.077, 3500;
2.879,3.036, 3249;3.421,3.793, 2149;3.953,3.036, 2849;3.402,2.219, 2849;
3.000,3.047, 3010;3.430,3.639, 3010;3.822,3.012, 2310;3.411,2.385, 3410;
3.103,3.012, 3345;3.430,3.462, 3845;3.710,3.036, 2645;3.402,2.562, 2745;
3.224,3.047, 3229;3.411,3.260, 3329;3.542,3.024, 3429;3.393,2.763, 3529];
x=A(:,1);y=A(:,2);z=A(:,3);
scatter(x,y,5,z)%散点图,5是点的大小
figure %打开显示图的界面
[X,Y,Z]=griddata(x,y,z,linspace(1.486,4.271)',linspace(1.604,4.276),'v4');%插值
pcolor(X,Y,Z);shading interp%伪彩色图
figure;contourf(X,Y,Z) %等高线图
figure;mesh(X,Y,Z)

A=[1.109,1.059,0.4874;2.021,0.841,0.5643;2.870,0.359,0.4628;4.039,0.196,0.4411;
4.505,3.012,0.4845;3.402,1.604,0.7857;3.570,3.565,0.7071;2.150,4.870,0.4284;
1.794,3.059,1.0000;2.121,3.615,0.8544;2.570,3.473,1.0000;3.421,4.160,0.5447;
4.271,3.036,0.5643;3.411,1.876,0.8771;2.561,2.562,1.0000;2.179,2.420,1.0000;
2.757,3.024,1.0000;3.439,3.970,0.6008;4.084,3.036,0.6325;3.402,2.077,0.9713;
2.879,3.036,1.0000;3.421,3.793,0.6667;3.953,3.036,0.6727;3.402,2.219,1.0000;
3.000,3.047,1.0000;3.430,3.639,0.7036;3.822,3.012,0.7180;3.411,4.215,0.5199;
1.103,4.612,0.3962;3.430,3.462,0.7857;3.710,3.036,0.7692;3.802,2.462,0.7670;
3.424,3.247,0.8771;3.511,3.060,0.8944;4.342,2.724,0.5522;3.803,2.903,0.7352];
x=A(:,1);y=A(:,2);z=A(:,3);
scatter(x,y,5,z)%散点图,5是点的大小
figure %打开显示图的界面
[X,Y,Z]=griddata(x,y,z,linspace(1.486,4.271)',linspace(1.604,4.276),'v4');%插值
pcolor(X,Y,Z);shading interp%伪彩色图
figure;contourf(X,Y,Z) %等高线图
figure;mesh(X,Y,Z)

matlab绘制温度场(尚未深入研究)

clear
echo on
d1=43;d2=7;dx=0.15;dy=0.1;xy=dx/dy;yx=dy/dx;
t=zeros(d1,d2);t1=ones(d1,d2);t0=zeros(d1,d2);
x=zeros(d1);y=zeros(d2);
x(1)=0;x(2)=dx/2;
for i=3:d1-1;
x(i)=x(i-1)+dx;
end
x(d1)=(d1-2)*dx;
y(1)=0;y(2)=dy/2;
for i=3:d2-1;
y(i)=y(i-1)+dy;
end
y(d2)=(d2-2)*dy;
t1=20*ones(d1,d2);
t=zeros(d1);
dt=0.1;ttt=30;
%nnn=tt/dt;
echo off
%for iii=1:nnn
% ttt=iii*dt;
tf=30;af=6.6;af=1/af;bta=-6.6;v=0.0625;
tin=100;tout=100;d=0.05;l=6;
for i=1:42;
if x(i)<v*ttt
t1(i,1)=tin+(300-tout)*x(i)/(v*ttt);
elseif x(i)>v*ttt+28*d
t1(i,1)=300-(300-tout)*(x(i)-v*ttt-28*d)/(l-28*d-v*ttt);
else
zz=-0.123*(x(i)-v*ttt)/d-3.52*exp(-0.123*(x(i)-v*ttt)/d);
t1(i,1)=10060*exp(zz);
end
end
for i=1:41;
for j=2:6
t1(i,j)=t1(i,1)-10*(j-1);
end
end
for iii=1:500;
t0=t1;
cd=1300;a=0.0003;
an=[1.69,-0.594,0.401,-0.168,0.027,-0.037,0.046,-0.05, 0.039,-0.012];
bn=[0, 0.333, 0.017,-0.131,0.054,0.0003,0.007,-0.012,0.026,0];
fncp=zeros(d1,d2);
for i=1:10;
fncp=an(i)*cos(i*pi/9*t1)+bn(i)*sin(i*pi/9*t1)+fncp;
end
fncp=fncp/4.1868;
b=1.3e-6;c=1.5e-9;
fnk=(a+b*t1+c*t1.^2)*360;
for i=2:42;
for j=2:6;
fnae(i,j)=2*yx*fnk(i,j)*fnk(i+1,j)/(fnk(i,j)+fnk(i+1,j));
fnaw(i,j)=2*yx*fnk(i,j)*fnk(i-1,j)/(fnk(i,j)+fnk(i-1,j));
fnan(i,j)=2*xy*fnk(i,j)*fnk(i,j+1)/(fnk(i,j)+fnk(i,j+1));
fnas(i,j)=2*xy*fnk(i,j)*fnk(i,j-1)/(fnk(i,j)+fnk(i,j-1));
end
end
fnap0=cd*fncp*dx*dy/dt;
fnbb=fnap0.*t1;
%for i=2:41
% for j=1:5
% t1(i,j)=t(i,1)-10*(j-1);
% end
% end
% t0=t1;
kk=af+0.5*dx/fnk(1,2);
bb=fnbb(2,2)+tf*dy/kk;
ap=fnae(2,2)+fnan(2,2)+fnap0(2,2)+dy/kk+2*xy*fnk(2,1)-bta*dx*dy;
fff=fnae(2,2)*t1(3,2)+fnan(2,2)*t1(2,3)+2*xy*fnk(2,1)*t1(2,1)+bb;
t1(2,2)=fff/ap;
for j=3:5;
kk=af+0.5*dx/fnk(1,j);
bb=fnbb(2,j)+tf*dy/kk;
ap=fnae(2,j)+fnas(2,j)+fnan(2,j)+fnap0(2,j)+dy/kk-bta*dx*dy;
fff=fnae(2,j)*t1(3,j)+fnan(2,j)*t1(2,j+1)+fnas(2,j)*t1(2,j-1)+bb;
t1(2,j)=fff/ap;
end
kk1=af+0.5*dx/fnk(1,6);
kk2=af+0.5*dy/fnk(1,6);
bb=fnbb(2,6)+tf*dy/kk1+tf*dx/kk2;
ap=fnae(2,6)+fnas(2,6)+fnap0(2,6)+dy/kk1+dx/kk2-bta*dx*dy;
fff=fnae(2,6)*t1(3,6)+fnas(2,6)*t1(2,5)+bb;
t1(2,6)=fff/ap;
for i=3:40;
for j=2:6;
if j==2;
as=2*xy*fnk(i,1);
fff=fnae(i,2)*t1(i+1,2)+fnaw(i,2)*t1(i-1,2)+fnan(i,2)*t1(i,3)+as*t1(i,1)+fnbb(i,2);
ap=fnae(i,2)+fnaw(i,2)+fnan(i,2)+as+fnap0(i,2)+dy/kk-bta*dx*dy;
t1(i,2)=fff/ap;
elseif j==6
kk=af+0.5*dy/fnk(i,6);
bb=fnbb(i,6)+tf*dx/kk;
fff=fnae(i,6)*t1(i+1,6)+fnaw(i,6)*t1(i-1,6)+fnas(i,6)*t1(i,4)+bb;
ap=fnae(i,6)+fnaw(i,6)+fnas(i,6)+fnap0(i,6)+dx/kk-bta*dx*dy;
t1(i,5)=fff/ap;
else
fff=fnae(i,j)*t1(i+1,j)+fnaw(i,j)*t1(i-1,j)+fnan(i,j)*t1(i,j+1)+fnas(i,j)*t1(i,j-1)+fnbb(i,j);
ap=fnae(i,j)+fnaw(i,j)+fnan(i,j)+fnas(i,j)+fnap0(i,j)-bta*dx*dy;
t1(i,j)=fff/ap;
end
end
end
for j=3:5;
kk=af+0.5*dx/fnk(42,j);
bb=fnbb(41,j)+tf*dy/kk;
ap=dy/kk+fnaw(41,j)+fnan(41,j)+fnas(41,j)+fnap0(41,j)-bta*dx*dy;
fff=fnaw(41,j)*t1(40,j)+fnan(41,j)*t1(41,j+1)+fnas(41,j)*t1(41,j-1)+bb;
t1(41,j)=fff/ap;
end
kk1=af+0.5*dx/fnk(42,6);
kk2=af+0.5*dy/fnk(41,7);
bb=fnbb(41,6)+tf*dy/kk1+tf*dx/kk2;
ap=fnaw(41,6)+fnas(41,6)+fnap0(41,6)+dy/kk1+dx/kk2-bta*dx*dy;
fff=fnaw(41,6)*t1(40,6)+fnas(41,6)*t1(41,5)+bb;
t1(41,6)=fff/ap;
kk=af+0.5*dx/fnk(42,2);
bb=fnbb(41,2)+tf*dy/kk;
ap=fnaw(41,2)+fnan(41,2)+fnap0(41,2)+dy/kk+2*xy*fnk(41,1)-bta*dx*dy;
fff=fnaw(41,2)*t1(40,2)+fnan(41,2)*t1(41,3)+2*xy*fnk(41,1)*t1(41,1)+bb;
t1(41,2)=fff/ap;

if abs(abs(t1)-abs(t0))<20;
break;
%else
%t1=t0+0.5*(t1-t0);
end
end
%if rem(iii,10)==0
figure
subplot(3,2,1)
t111=t1';
%mesh(t111);
%view(2);
pcolor(t111);
subplot(3,2,5)
ccc=contour(t111,[200,400,600,800,1000]);
%ccc=contour(t111,5);
clabel(ccc);
%end
% 将平面网格映射到三维网格面上
x1=1:7;
y1=1:43;
[X1,Y1]=meshgrid(y1,x1);
%绘制加密前三维网格
subplot(3,2,3);
mesh(X1,Y1,t111);
% 用二元差值公式将网格加密并绘出温度场的三维图
x2=1:.25:7;
y2=1:.25:43;
[X2,Y2]=meshgrid(y2,x2);
Z2=griddata(X1,Y1,t111,X2,Y2);
t111=round(t111);
subplot(3,2,4);
mesh(X2,Y2,Z2)
% 绘制加密后的平面网格
subplot(3,2,2);
pcolor(Z2)
% 绘制加密后的温度场等高线图
subplot(3,2,6)
CCC=contour(Z2,[100,200,300,400,500,600,700,800,900,1000]);
clabel(CCC);
% end
%下面输入平板中某点坐标以查看此点温度
a1=input('请输入你要查看温度的点的X轴坐标值（范围为1到169）');
b1=input('请输入你要查看温度的点的Y轴坐标值（范围为1到25）');
disp('输入坐标值为'),a1,b1
a2=round((a1+3)/4);
b2=round((b1+3)/4);
disp('计算出网格加密前坐标点的温度为')
Q1=(t111(b2,a2))
display('计算出网格加密后坐标点的温度为')
Q2= Z2(b1,a1)

二维曲线（非线性）拟合步骤

1. 建立函数
function F = myfun(x,xdata)
%F = x(1)*xdata.^2 + x(2)*sin(xdata) + x(3)*xdata.^3;
F = x(1)+x(2)*xdata.^2 + x(3)*sin(xdata) + x(4)*xdata.^3; % 这个函数的可以是任意的,这个要%根据实际情况自己来选;

2.然后给出数据xdata和ydata
xdata = -1.0:0.5:2.0;
ydata = [-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];
x0 = [10, 10, 10,1]; %初始估计值,可以随意给,有了这个初值后matlab才能进行迭代, 这是matlab的迭代的初始条件, 这个初始条件与结果无关; % [ ]中的初值的数量要和myfun中x(1) x(2) x(3) x(4)的数量一致才行
%拟合函数的表达式:
[x,resnorm] = lsqcurvefit(@myfun,x0,xdata,ydata) %回车后;
x =
0.5491 -2.9977 -0.0000 1.9991
x结果即为myfun 中对应的x(1) x(2) x(3) x(4)的值；这样就得到了所求的多项式的系数, resnorm是在x处残差的平方和, 在workpsace中可以查到, resnorm=sum ((fun(x,xdata)-ydata).^2)，
函数建立成功后,在matlab输入的程序如下;
clear;%可以有,也可以没有,最好有.
xdata = -1.0:0.5:2.0;
ydata = [-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];
x0 = [10, 10, 10, 1]; %初始估计值
[x,resnorm] = lsqcurvefit(@myfun,x0,xdata,ydata)
怎么样,很简单吧?
如果不想编写函数的话也可以将最后一句话改为:
[x,resnorm] = lsqcurvefit(@(x,xdata)x(1)+x(2)*xdata.^2 + x(3)*sin(xdata) + x(4)*xdata.^3,x0,xdata,ydata);这样就不用编写函数了;

三维曲线（非线性）拟合步骤

1 设定目标函数. （M函数书写）% 可以是任意的
例如：

function f=mydata(a,data) %y的值目标函数值 或者是第三维的,a=[a(1) ，a(2)] 列向量
x=data(1,:); %data 是一2维数组，x=x1
y=data(2,:); %data 是一2维数组，x=x2
f=a(1)*x+a(2)*x.*y; %这里的a(1)， a(2)为目标函数的系数值。 f的值相当于ydata的值

2 然后给出数据xdata和ydata的数据和拟合函数lsqcurvefit

例如：
x1=[1.0500 1.0520 1.0530 1.0900 1.0990 1.1020 1.1240 1.1420...
1.1490 1.0500 1.0520 1.0530 1.0900 1.0990 1.1020 1.1240 1.1420 1.1490];
x2=[3.8500 1.6500 2.7500 5.5000 7.7000 3.3000 4.9500 8.2500 11.5500...
1.6500 2.7500 3.8500 7.7000 3.3000 5.5000 8.2500 11.5500 4.9500];
ydata=[56.2000 62.8000 62.2000 40.8000 61.4000 57.5000 44.5000 54.8000...
53.9000 64.2000 62.9000 64.1000 63.0000 62.2000 64.2000 63.6000...
52.5000 62.0000];
data=[x1;x2]; %类似于将x1 x2整合成一个2维数组。
a0= [-0.0014,0.07];
option=optimset('MaxFunEvals',5000);
format long;
[a,resnorm]=lsqcurvefit(@mydata,a0,data,ydata,[],[],option);
yy=mydata(a,data);
result=[ydata' yy' (yy-ydata)']
将数据按列的形式排列:第一列是yadata; 第二列是yadata(根据拟合函数得到的计算值);第三列是上述二者的插值; ' 这个符号实现了数据按列进行排列的功能;这个符号后面有空格
result =
56.200000000000003 57.919539361940053 1.719539361940051
62.799999999999997 60.487127691642989 -2.312872308357008
62.200000000000003 59.314824360714432 -2.885175639285571
40.799999999999997 58.216478553229628 17.416478553229631
61.399999999999999 56.130116906144998 -5.269883093855000
57.500000000000000 61.431449491527204 3.931449491527204
44.500000000000000 60.688766221557501 16.188766221557501
54.799999999999997 57.659418794420404 2.859418794420407
53.899999999999999 53.987090284681393 0.087090284681395
64.200000000000003 60.372133152305260 -3.827866847694743
62.899999999999999 59.258494992850515 -3.641505007149483
64.099999999999994 58.085023760117025 -6.014976239882969
63.000000000000000 55.670452618469561 -7.329547381530439
62.200000000000003 61.264213240642817 -0.935786759357185
64.200000000000003 58.857393913448675 -5.342606086551328
63.600000000000001 56.750601335313952 -6.849398664686049
52.500000000000000 53.658187210710317 1.158187210710317
62.000000000000000 62.038605327908876 0.038605327908876
% a的值为拟合的目标函数的参数值 利用lsqcurvefit进行拟合的 它完整的语法形式是：
% [x,resnorm,residual,exitflag,output,lambda,jacobian] =lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)

三维曲线的画法
三维空间