第2章 matlab.docx
《第2章 matlab.docx》由会员分享，可在线阅读，更多相关《第2章 matlab.docx（17页珍藏版）》请在冰豆网上搜索。
第2章matlab
第二章MATLAB绘图与插值拟合
§2.2二维作图
1.数值图
例2.2.3:
t=0:
0.1:
2*pi;plotyy（t,sin（t）,t,0.01*cos（t））
例2.2.4
%li2_2_4.m
n=1:
100;
y1=（1+1./n）.^n;y=（1+1./n）.^（n+1）;
plot（n,y1,'.',n,y,'.r'）
例2.2.5
function[x1,y1]=li2_2_5（n）
%以li2_2_5.m文件名存盘
t1=1;x
（1）=1;t2=2;y
（1）=2;
fork=1:
n
x（k+1）=sqrt（x（k）*y（k））;y（k+1）=（x（k）+y（k））/2;
end
x1=x（k+1）;
y1=y（k+1）;
t=1:
n+1;
plot（t,x,'-r',t,y,'.b'）
例2.2.6
%li2_2_6.m
t=0:
0.1:
2*pi;y1=sin（t）;y2=cos（t）;y3=sin（t）*cos（t）;
plot（t,y1,'-.g',t,y2,':
r',t,y3,'xb'）
例2.2.7
%li2_2_7.m
x=-pi:
pi/10:
pi;y=tan（sin（x））-sin（tan（x））;
plot（x,y,'--ro','LineWidth',2,'MarkerEdgeColor','k',…
'MarkerFaceColor','g','MarkerSize',10）
例2.2.8theta=0:
0.1:
8*pi;polar（theta,cos（5*theta）+1/3）
例2.2.9
%li2_2_9.m
x=-2:
0.2:
2;y=sin（x）;
subplot（2,2,1）,plot（x,y）
subplot（2,2,3）,hist（y,5）
subplot（2,2,4）,stem（x,y）;
subplot（2,2,2）,rose（x,y）;%stem（）和rose也是绘制二维曲线的函数。
Subplot（4,4,11）,fill（x,y,'r'）;
Subplot（4,4,12）,feather（x,y）;
Subplot（4,4,15）,plot（x,y）
Subplot（4,4,16）,stairs（x,y）;
例2.2.10
%li2_2_10.m
subplot（2,2,1）,
ezplot（'5*x^4+2*x^2+1'）
subplot（2,2,2）,ezplot（'sin（x）*exp（-x^2）'）
subplot（2,2,3）,
ezplot（'cos（x）/cos（2*x）'）
subplot（2,2,4）,
ezplot（'1+36*x/（x+3）^2'）
§2.3三维图形的绘制
例2.3.1
%li2_3_1.m
t=0:
pi/50:
10*pi
x=sin（t）;y=cos（t）;z=t;h=plot3（x,y,z,'g-'）
set（h,'LineWidth',4*get（h,'LineWidth'））
例2.3.2
%li2_3_2.m
[x,y]=meshgrid（-3:
0.1:
3,-3:
0.1:
3）;%生成平面网格
z=x.^2-y.^2;
mesh（x,y,z）
例2.3.3
[x,y]=meshgrid（-3:
0.1:
3,-2:
0.1:
2）;%li2_3_3.m
z=（x.^2-2*x）.*exp（-x.^2-y.^2-x.*y）;
axis（[-3,3,-2,2,-0.7,1.5]）%设定坐标系范围
mesh（x,y,z）
§2.4车灯光源投影区域的绘制（CUMCM2002A案例）
%li2_4_1.m
p=0.03;x=25.0216;
fory1=-0.002:
0.0004:
0.002
y0=（-0.036:
0.001:
0.036）'*ones（1,73）;
z0=ones（73,1）*（-0.036:
0.001:
0.036）;
x0=（y0.^2+z0.^2）/（2*p）;
xn=（p^3+4*x0*2*p.*x0+p*（-4*y1*y0+3*2*p*x0））./（2*（p^2+2*p*x0））;
yn=（2*p*x0.*y0+p^2*（-y1+y0）+y1*（y0.^2-z0.^2））./（p^2+2*p*x0）;
zn=（p^2+2*p*x0+2*y1*y0）.*z0./（p^2+2*p*x0）;
y=y0+（yn-y0）.*（x-x0）./（xn-x0）;
z=z0+（zn-z0）.*（x-x0）./（xn-x0）;
plot（y,z,'b.'）
xlabel（'y'）
ylabel（'z'）
holdon
end
§2.5动画的绘制
例2.5.1
%li2_5_1.m
clear
a=[-8/300;0-1010;028-1];
y=[35-10-7]';
h=0.01;
p=plot3（y
（1）,y
（2）,y（3）,'.','EraseMode','none','MarkerSize',5）;
%将擦除模式设置为none
axis（[050-2525-2525]）
holdon
fori=1:
4000
a（1,3）=y
（2）;
a（3,1）=-y
（2）;
ydot=a*y;
y=y+h*ydot;
set（p,'XData',y
（1）,'YData',y
（2）,'ZData',y（3））%设定图形目标的性质
drawnow%填充未完成的图形事件
i=i+1;
end
例2.5.2
%li2_5_2.m
clear
axis（[0,2*pi,-1,1]）
x=0:
0.01:
2*pi;
plot（x,0）
holdoff
holdon
forj=1:
30
t1=（j-1）*2*pi/30;
t2=j*2*pi/30;
t=t1:
0.01:
t2;
plot（t,0,'.r'）
plot（t,sin（t）,'.'）
end
例2.5.3
%li2_5_3.m
clear
axisequal%创建一个用于显示该电影动画的坐标轴
m=moviein（16）;
set（gca,'NextPlot','replacechildren'）
forj=1:
16
plot（fft（eye（j+16）））
m（:
j）=getframe;
end
movie（m,30）
例2.5.4
%li2_5_4.m
clearall
axis（[0,2*pi,-1,1]）
holdoff
x=0:
0.01:
2*pi;
plot（x,0）
m=moviein（16）;
holdon
fori=1:
50
forj=1:
i
t1=（j-1）*2*pi/50;
t2=j*2*pi/50;
t=t1:
0.01:
t2;
axis（[0,2*pi,-1,1]）%限制动画的坐标显示大小
plot（t,0,'.r'）
plot（t,sin（t）,'.'）
end
m（:
i）=getframe;
end
movie（m,30）
例2.5.5
functionx1=li2_5_5（a）
%li2_5_5.m作函数
极限的
语言的电影动画程序
%a取0.1~0.01之间
holdoff
axis（[0,70,0,1.8]）
x=0.01:
0.01:
70;
y=（x-sin（x））./（x+sin（x））;
plot（x,y,'-g'）;
holdon
plot（x,1+a,'-b'）;
plot（x,1-a,'-b'）;
m=moviein（16）;
fork=70:
-0.1:
1
ifabs（（（k-sin（k））/（k+sin（k）））-1）>a
x1=k;
break;
end
end
n2=floor（x1）+1;
forj=1:
n2
ifj==1
t1=0.01;
t2=1;
else
t1=j-1;
t2=j;
end
t=t1:
0.01:
t2;
axis（[0,70,0,1.8]）%设定坐标轴显示尺寸
y=（t-sin（t））./（t+sin（t））;
plot（t,0,'.r'）;
plot（t,y,'-r'）;
m（:
j）=getframe;
end
t=0:
0.01:
1+a;
n=length（t）;
n1=n2*ones（1,n）;
plot（n1,t,'.g'）
m（:
n2+1）=getframe;
movie（m,30）
§2.6插值
例2.6.1
functiony=li2_6_1（x0,y0,x）
%li2_6_1.m拉格朗日插值函数
n=length（x0）;m=length（x）;
fori=1:
m
z=x（i）;
s=0;
fork=1:
n
p=1;
forj=1:
n
ifj~=k
p=p*（z－x0（j））/（x0（k）－x0（j））；
end
end
s=p*y0（k）＋s;
end
y（i）=s;
end
%range.m
n=11;m=21;
x0=
-5:
10/（n-1）:
5;
y0=1./（1+x0.^2）;
x=-5:
10/（m-1）:
5;
y=1./（1+x.^2）;
y1=shiyan7_1（x0,y0,x）;
y2=interp1（x0,y0,x）;
y3=interp1（x0,y0,x,'spline'）;
s='xyy1y2y3'
[x（11:
21）'y（11:
21）'y1（11:
21）'y2（11:
21）'y3（11:
21）']
z=0*x；
plot（x,z,x,y,'k:
',x,y2,'r','LineWidth',1.5）
figure
（2）,
plot（x,z,x,y,'k:
',x,y3,'r','LineWidth',1.5）
例2.6.2
%lt2_6_2a.m
x=1:
5;
y=1:
3;
temps=[8281808284;7963616581;8484828586];
mesh（x,y,temps）
%lt2_6_2b.m
x=1:
5;
y=1:
3;
temps=[8281808284;7963616581;8484828586];
cx=1:
0.05:
5;
cy=1:
0.05:
3;
cz=interp2（x,y,temps,cx',cy,'cubic'）;
mesh（cx,cy,cz）
%lt2_6_2c.m
x=[129140103.588185.5195105.5157.5107.57781162162117.5];
y=[7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5];
plot（x,y,'+'）
作出海底地貌图
编程lt2_6_2d.m如下：
%lt2_6_2d.m
x=[129140103.588185.5195105.5157.5107.57781162162117.5];
y=[7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5];
z=[-4-8-6-8-6-8-8-9-9-8-8-9-4-9];
cx=75:
0.5:
200;
cy=-70:
0.5:
150;
cz=griddata（x,y,z,cx,cy','v4'）;
meshz（cx,cy,cz）,rotate3d%作海底地貌图
xlabel（'X'）,ylabel（'Y'）,zlabel（'Z'）
figure
（2）,contour（cx,cy,cz,[-5-5]）;grid%作深度为5的海底等值线图
xlabel（'X'）,ylabel（'Y'）
[i1,j1]=find（cz<-5）;
fork=1:
length（i1）
cz（i1（k）,j1（k））=-5;
end
figure（3）,meshc（cx,cy,cz）,rotate3d%作水深低于5英尺的部分海底曲面图
例2.6.3
%lt2_6_3.m
[x,y]=meshgrid（1:
10）;%构造测量网络
h=[0,0.02,-0.12,0,-2.09,0,-0.58,-0.08,0,0;...
0.02,0,0,-2.38,0,-4.96,0,0,0,-0.1;...
0,0.1,1,0,-3.04,0,-0.53,0,0.1,0;...
0,0,0,3.52,0,0,0,0,0,0;...
-0.43,-1.98,0,0,0,0.77,0,2.17,0,0;...
0,0,-2.29,0,0.69,0,2.59,0,0.3,0;...
-0.09,-0.31,0,0,0,4.27,0,0,0,-0.01;...
0,0,0,5.13,7.4,0,1.89,0,0.04,0;...
0.1,0,0.58,0,0,1.75,0,-0.11,0,0;...
0,-0.01,0,0,0.3,0,0,0,0,0.01];%测量数据点[xi,yi]=meshgrid（1:
0.1:
10）;%构造插值网络
hi=interp2（x,y,h,xi,yi,'spline'）;%二维插值命令
surf（hi）;%画出地貌图
xlabel（'x'）,ylabel（'y'）,zlabel（'h'）;%添加坐标名称
%li2_7_1.m重金属污染地形图
A=xlsread（'F:
\dimao.xls'）;读取F盘上地貌数据文件dimao.xls
x=A（:
1）;
y=A（:
2）;
z=A（:
3）;
scatter（x,y,5,z）%散点图
figure
[X,Y,Z]=griddata（x,y,z,linspace（min（x）,max（x）,200）',linspace（min（y）,max（y）,200）,'v4'）;%插值
figure,surf（X,Y,Z,'EdgeColor','none'）%三维曲面
title（‘重金属污染地形图’）
grifoff
%li2_7_2.m重金属浓度随地形分布图
clearall
A=xlsread（'F:
\dimao.xls'）;
B=xlsread（'F:
\zjs.xls'）;
x=A（:
1）;
y=A（:
2）;
z=B（:
1）;%As数据，第i个重金属数据对应B矩阵中第i列
scatter（x,y,5,z）%散点图
figure
[X,Y,Z]=griddata（x,y,z,linspace（min（x）,max（x）,200）',linspace（min（y）,max（y）,200）,'v4'）;%插值
figure,surf（X,Y,Z,'EdgeColor','none'）%三维曲面
title（'As浓度随地域分布图'）%As随对应重金属所在列做改动
grifoff
§2.8拟合
例2.8.1编写如下程序
x=[0.5,1.0,1.5,2.0,2.5,3.0]%给出数据点的x值
y=[1.75，2.45,3.81,4.80,7.00,8.60]；%给出数据点的y值
p=polyfit（x,y,2）；%求出2阶拟合多项式的系数
poly2str（p,＇x＇）%求出2次拟和多项式表达式
＞＞x=-2*pi:
0.2:
2*pi；y=1/5*cosh（x./3）；%将自变量和函数离散化
＞＞p=polyfit（x,y,2）%求出拟合多项式的系数
P=
0.0151-0.0003-.1836
＞＞pk=poly2str（p,＇x＇）%求出拟合多项式表达式
pk=
0.015071x^2-0.00029557x+0.18365
＞＞ezpolt1/5*cosh（x/3），hold%作悬链线函数曲线并保留
＞＞ezpoto.o15071*x^2-0.00029557*x+0.18365%作拟合多项式函数曲线
＞＞grid,legend（＇悬链线函数线＇，＇拟合悬链线的二次多项式曲线＇）
例2.8.3
function[a,b]=li2_8_3（c1,c2）;
%存盘文件名li2_8_3.m
%做非线性曲线拟合，参数传递时取c1=2;c2=1;
x=1:
5;
y=c1*exp（c2*x）+20*normrnd（0,1,1,5）;%在指数曲线上加上随机误差，产生5个散点
[a,b]=lsqcurvefit（@nh,[c1,c2],x,y）;%非线性曲线拟合指令调用nh.m函数
x1=1:
0.01:
5;
y2=c1*exp（c2*x1）;
y3=a
（1）*exp（a
（2）*x1）;
plot（x1,y2,'b-',x,y,'ro',x1,y3,'g-'）;%作散点图、原指数（
）图、拟合图
End
functionf=nh（a,x）%x数据通过主函数传递，a是由主函数传递过来的c1,c2迭代初值
%子函数程序，存盘文件名nh.m
f=a
（1）*exp（a
（2）.*x）;
end
例2.8.4
li2_8_4a.m.
functionf=fun1（a,x）
%li2_8_4a.m.
f=exp（a
（1）*x+a
（2））;%用最小乘二拟合求上述函数中待定常数，
%以及检验拟合效果的图形绘制程序
t=1790:
10:
1990;
x=[3.95.37.29.612.917.123.231.438.650.262.976…
92106.5123.2131.7150.7179.3204226.5251.4];
plot（t,x,’*’,t,x）;
a0=[0.001,1];
a=lsqcurvefit（‘fun1’,a0,t,x）
ti=1790:
10:
2020;
xi=fun1（a,ti）;
holdon
plot（ti,xi）;
t1=2010;
x1=fun1（a,t1）
holdoff
阻滞增长模型（logistic模型）.
在MATLAB命令窗口键入：
Dsolve（’Dx=r*x*（1-x/c）’,’x（1790）=3.9’）
ans=
c/（1+1/39*exp（-r*t）*exp（r）^1790*（10*c-39））
定义函数（2-8-3）的函数M文件.
functionf=fun3（a,t）
f=a
（1）./（1+（a
（1）/3.9-1）*exp（-（t-1790）*a
（2）））;
%li2_8_4b.m.
x=1790:
10:
1990;
y=[3.95.37.29.612.917.123.231.438.650.262.976...
92106.5123.2131.7150.7179.3204206.5251.4];
plot（x,y,’*’,x,y）;
a0=[0.001,1];
a=lsqcurvefit（‘fun3’,a0,x,y）
xi=1790:
5:
2020;yi=fun3（a,xi）;
holdon
plot（xi,yi）;
x1=2010;
y1=fun3（a,x1）
先编写一个函数M文件，以定义优化问题（2-8-4）式中的目标函数。
functionf=fun5（a）
n=16;w=2;
x=1790:
10:
1990;x1=x（1:
n）;x2=x（n+1:
21）;
y=[3.95.37.29.612.917.123.231.438.650.262.976...
92106.5123.2131.7150.7179.3204226.5251.4];
y1=y（1:
n）;y2=y（n+1:
21）;
f=[fun3（a,x1）,w*fun3（a,x2）]-[y1,w*y2];
主程序：
%li2_8_4c.m.
t=1790:
10:
1990;
x=[3.95.37.29.612.917.123.231.438.650.262.976...
92106.5123.2131.7150.7179.3204226.5251.4];
plot（t,x,'*',t,x）;
a0=[300,0.03];
a=lsqnonlin（'fun5',a0）
ti=1790:
5:
1990;
xi=fun3（a,ti）;
holdon;
plot（ti,xi,t,x,'*'）;
x1=fun3（a,2010）
holdoff
例2.8.5
%Li2_8_5.m
clf,clear,
n=10;m=3;x=1:
1:
12;
y=[3.13.86.912.716.820.524.525.922.016.110.75.4];
z=（pi/6）*x;plot（z（1:
12）,y（1:
12）,'o'）;holdon
k=1:
m;%计算数据矩阵.
fori=1:
n
X（i,:
）=[1cos（k*z（i））sin（k*z（i））];
end
c=inv（X'*X）*X'*y（1:
n）',%求解.
z1=linspace（0,2*pi,30）;
s=[];%开始计算F-级数的部分和.
fori=1:
30;
sd=[1cos（k*z1（i））sin（k*z1（i））]';%注意k是向量.
sd=c.*sd;s=[s,sum（sd）];
end
plot（z1,s,'r-'）,holdon,xlabel（'月份'）,ylabel（'平均气温'）
§2.9黄河小浪底调水调沙问题
clc,clearall%Li2_9_1a.m
Loaddata3.txt%把表中5.8中的日期和时间数据行删除，余下的数据保存在纯文本文件
Liu=data3（[1:
3,:
]）;liu=liu’;liu=liu（:
）;%提出水流量并按照顺序变成列向量
sha=data3（[2:
4,:
]）;sha=sha’;sha=sha（:
）;%提出含沙量并按照顺序变成列向量
y=sha.*liu;y=y’;%计算排沙量，并变成行向量
i=1:
24；
t=（12*i-4）*3600;
t1=t
（1）;t2=t（end）;
pp=csape（t,p）;%进行三次样条差值
xsh=pp.coefs%求得差值多项式的系数矩阵，每一
升级会员 