数学实验与matlab 程序解析资料全编.docx
《数学实验与matlab 程序解析资料全编.docx》由会员分享，可在线阅读，更多相关《数学实验与matlab 程序解析资料全编.docx（85页珍藏版）》请在冰豆网上搜索。
数学实验与matlab程序解析资料全编
《数学实验与Matlab》程序
将程序按书中的顺序排列，本书程序均通过了调式。
可直接拷贝到命令窗口运行或编制M文件运行。
如出现问题，可能是中英文标点的缘故（排版有可能使用了中文标点），请将中文标点换为英文标点试试。
实验１．矩阵运算与Matlab命令
1.1知识要点与背景：
知识要点和背景：
日常矩阵及其运算
【A=[423;132;133;322],%表1-1、表1-2的数据分别写成矩阵形式
B=[35206045;10155040;20124520]】
【C=A*B%矩阵乘法，求各订单所对应的原材料和劳动力。
】
【whos%查看Matlab工作空间中变量及其规模】
1.2实验与观察：
矩阵和Matlab语言
1.2.1向量的生成和运算
【x=linspace（0,4*pi,100）;%将[0,4π]区间100等分，产生了一个100维向量
y=sin（x）;%计算函数值，产生了一个与x同维的100维函数向量y
y1=sin（x）.^2;%计算函数向量，注意元素群运算
y2=exp（-x）.*sin（x）;
%以x为横坐标，y为纵坐标画函数的图用不同的线型将函数曲线绘制在一个图上
plot（x,y,'-',x,y1,'-',x,y2,'.-'）】
1.向量的创建
◆直接输入向量。
【x1=[124],x2=[1,2,1],x3=x1'】
◆冒号创建向量。
【x1=3.4:
6.7
x2=3.4:
2:
6.7
x3=2.6:
-0.8:
0】
◆生成线性等分向量。
【x=linspace（0,1,5）】
2.向量的运算
【y=sin（x）】
【y1=sin（x）.^2;y2=exp（-x）.*sin（x）;】
1.2.2.矩阵创建和运算
1.创建矩阵
（1）数值矩阵的创建
◆直接输入法创建简单矩阵。
【A=[1234;5678;9101112]】
【B=[-1.3,sqrt（3）;（1+2）*4/5,sin（5）;exp
（2）,6]】
（2）符号矩阵的创建
◆
【symsa11a12a13a14a21a22a23a24a31a32a33a34…
b11b12b13b14b21b22b23b24b31b32b33b34
A1=[a11a12a13a14;a21a22a23a24;a31a32a33a34],
B1=[b11b12b13b14;b21b22b23b24;b31b32b33b34]】
2.矩阵的运算
【C=A1+B1,D=A1-B1】
【symsc
cA=c*A1】
【C=A1*B1】
{?
?
?
Errorusing==>sym/mtimes,Innermatrixdimensionsmustagree.}
【A2=A1（:
1:
3）,B1】
【G=A2*B1】
【g11=A2（1,:
）*B1（:
1）】
【A,A_trans=A'】
【H=[123;210;123],K=[123;210;231]
h_det=det（H）,k_det=det（K）,H_inv=inv（H）,K_inv=K^-1】
【A=[301;110;014];
B=inv（A-2*eye（3））*A,B=（A-2*eye（3））\A】
3.分块矩阵：
矩阵的裁剪、分割、修改与抽取
（1）
【A=[10112;01-123;30510;23121],vr=[1,3];vc=[1,3];
A1=A（vr,vc）%取出A的1、3行和1、3列的交叉处元素构成新矩阵A1】
◆将上面的矩阵A分为四块，并把它们赋值到矩阵B中，观察运行后的结果。
【A11=A（1:
2,1:
2）,A12=A（1:
2,3:
5）,A21=A（3:
4,1:
2）,A22=A（3:
4,3:
5）
B=[A11A12;A21A22]】
A=
205
421
0-12
B=
124-1
5310
-1023
C=
-341813
131420-1
-7-336
（2）矩阵的修改和提取
◆【A=[10112;01-123;30510;23121]
A（1,:
）=[00000];A】
◆观察：
【B（:
[2,4]）=[]%删除矩阵B的第2、4列】
（3）矩阵元素的抽取
4.生成特殊矩阵
。
◆
【y1=rand（1,5）,y2=rand（1,5）,
rand（'seed',3）,x1=rand（1,5）,rand（'seed',3）,x2=rand（1,5）】
5.常用矩阵函数
6.数据的简单分析
◆
【rand（'seed',1）;A=rand（3,6）,
Asort=sort（A）,Amax=max（A）,Asum=sum（A）】
1.2.3Matlab工作环境和编程
2.Matlab的基本设计
1.3应用、思考与练习
1.3.1关系矩阵
1.3.2投入产出
1.3.3循环比赛的名次
【A=[0110;0011;0001;1000],
e=ones（4,1）;c=A*e;
s=c'】
★画矩阵结构图的gplot指令。
◆（3）
【clf,A=[0110;0011;0001;1000];xy=[01;00;-1–0.5;1–0.5];
graphy_plot（A,xy,1,0.5）,%gplot（A,xy）】
1.3.4参考程序
graphy_plot.m
【functiony=graphy_plot（A,xy,l,p）
%画矩阵的有向结构图。
A为邻接矩阵，xy为顶点坐标，l控制参数，l=0，画无向图；
%l~=0,画有向图。
p为控制箭头大小的参数。
a=-max（abs（xy（:
1）））*1.1;b=max（abs（xy（:
1）））*1.1;
c=-max（abs（xy（:
2）））*1.1;d=max（abs（xy（:
2）））*1.1;
ifl=0
gplot（A,xy）,axis（[abcd]）,holdon,
elseifl~=0
U=[];V=[];X=[];Y=[];
n=length（A（:
1））;
fori=1:
n
k=find（A（i,:
）~=0）;
m=length（k）;
if（m~=0）
forj=1:
m
u
（1）=（xy（k（j）,1）-xy（i,1））;v
（1）=（xy（k（j）,2）-xy（i,2））;
u
（2）=eps;v
（2）=eps;U=[u;U];V=[v;V];
X=[[xy（i,1）xy（k（j）,1）];X];Y=[[xy（i,2）xy（k（j）,2）];Y];
end
text（xy（i,1）,xy（i,2）,['\bullet\leftarrow\fontsize{16}\it{V}',…
um2str（i）]）;holdon,
end
end
gplot（A,xy）,axis（[abcd]）,holdon,
h=quiver（X,Y,U,V,p）;set（h,'color','red'）;holdon,
plot（xy（:
1）,xy（:
2）,'k.','markersize',12）,holdon,
end,holdoff】
实验２．函数的可视化与Matlab作
2.1实验与观察：
函数的可视化
2.1.1Matlab二维绘图命令
1.周期函数与线性p－周期函数
◆观察：
【clf,x=linspace（0,8*pi,100）;
F=inline（'sin（x+cos（x+sin（x）））'）;
y1=sin（x+cos（x+sin（x）））;y2=0.2*x+sin（x+cos（x+sin（x）））;
plot（x,y1,'k:
',x,y2,'k-'）legend（'sin（x+cos（x+sin（x））','0.2x+sin（x+cos（x+sin（x）））',2）】
2.plot指令：
绘制直角坐标的二维曲线
3.图形的属性设置和屏幕控制
【h=plot（[0:
0.1:
2*pi],sin（[0:
0.1:
2*pi]））;gridon
set（h,'LineWidth',5,'color','red'）;set（gca,'GridLineStyle','-','fontsize',16）】
◆设置y坐标的刻度并加以说明，并改变字体的大小。
【h=plot（[0:
0.1:
2*pi],sin（[0:
0.1:
2*pi]））;gridon，
set（gca,'ytick',[-1-0.500.51]）,set（gca,'yticklabel','a|b|c|d|e'）,
set（gca,'fontsize',20）】
4.文字标注指令
【plot（x,y1,'b',x,y2,'k-'）,
set（gca,'fontsize',15,'fontname','timesNewRoman'）,%设置轴对象的字体为times
%NewRoman，字体的大小为15
title（'\it{Peroidandlinearperoidfunction}'）;%加标题，注意文字用单引号''加上
%斜杠'\'后可输入不同的设置，例如it{…}表示花括号里的文字为斜体；如果有多项设置，
%则可用\…\…\…连续输入。
xlabel（'xfrom0to8*piit{t}\'）;ylabel（'\it{y}'）;%说明坐标轴
text（x（49）,y1（50）-0.4,'\fontsize{15}\bullet\leftarrowTheperiodfunction{\itf（x）}'）;
%在坐标（x（49）,y1（50）-0.4）处作文字说明，各项设置用"\"隔开。
%\fontsize{15}\bullet\leftarrow的意义依次是：
\字体大小=15\画圆点\左箭头
text（x（14）,y2（50）+1,'\fontsize{15}Thelinearperiodfunction{\itg（x）}\rightarrow
\bullet'）%与上一语句类似，用右箭头】
图2.5文字标注
◆观察指令legend和num2str的用法：
在同一张图上画出
，这里
并进行适当的标注。
zxy2_2.m
【clf,t=0:
0.1:
3*pi;alpha=0:
0.1:
3*pi;
plot（t,sin（t）,'r-'）;holdon;plot（alpha,3*exp（-0.5*alpha）,'k:
'）;
set（gca,'fontsize',15,'fontname','timesNewRoman'）,
xlabel（'\it{t（deg）}'）;ylabel（'\it{magnitude}'）;
title（'\it{sinewaveand{\it{Ae}}^{-\alpha{\itt}}wave}'）;%注意\alpha的意义
text（6,sin（6）,'\fontsize{15}TheValue\it{sin（t）}at{\itt}=6\rightarrow\bullet','HorizontalAlignment','right'）,
%上面的语句是一整行，如果要写成两行，必须使用续行号…，例如要在“bullet',”
%后换行，需写“bullet',…”后才能换行。
%'HorizontalAlignment','right'表示箭头所指的曲线对象在文字的右边。
text（2,3*exp（-0.5*2）,['\fontsize{15}\bullet\leftarrowTheValueof\it{3e}^{-0.5\it{t}}=',num2str（3*exp（-0.5*2））,'at\it{t}=2']）;
%num2str的用法：
['string1',num2str,'string2']，注意方括号的使用。
%legend（'\itsin（t）','{\itAe}^{-\alphat}'）%请结合图形观察此命令的使用】
运行结果如图2.6所示。
5.图形窗口的创建和分割
◆观察：
【clf，b=2*pi;x=linspace（0,b,50）;
fork=1:
9
y=sin（k*x）;
subplot（3,3,k）,plot（x,y）,axis（[0,2*pi,-1,1]）
end】
x=linspace（0,4*pi,150）;
y1=cos（x）;
y=cos（cos（cos（cos（cos（cos（cos（x）+sin（x）））+sin（x）））+sin（x）））;
plot（x,y,x,y1,'r-'）,grid
2.1.2多元函数的可视化与空间解析几何（三维图形）
本节通过高等数学的几个例子观察Matlab的三维绘图功能和技巧。
1.绘制二元函数
◆观察：
绘制
的图象，作定义域的裁剪。
◆
（1）观察meshgrid指令的效果。
【a=-0.98;b=0.98;c=-1;d=1;n=10;
x=linspace（a,b,n）;y=linspace（c,d,n）;
[X,Y]=meshgrid（x,y）;
plot（X,Y,'+'）】
★三维绘图指令mesh、meshc、surf。
◆
（2）做函数的定义域裁剪，观察上述三维绘图指令的效果。
程序zxy2_4.m
【clear,clf,
a=-1;b=1;c=-15;d=15;n=20;eps1=0.01;
x=linspace（a,b,n）;y=linspace（c,d,n）;
[X,Y]=meshgrid（x,y）;
fori=1:
n%计算函数值z，并作定义域裁剪
forj=1:
n
if（1-X（i,j））z（i,j）=NaN;%作定义域裁剪，定义域以外的函数值为NaN
else
z（i,j）=1000*sqrt（1-X（i,j））^-1.*log（X（i,j）-Y（i,j））;
end
end
end
zz=-20*ones（1,n）;plot3（x,x,zz）,gridoff,holdon%画定义域的边界线
mesh（X,Y,z）%绘图，读者可用meshz,surf,meshc在此替换之
xlabel（'x'）,ylabel（'y'）,zlabel（'z'）,boxon%把三维图形封闭在箱体里】
◆运行了zxy2_4.m以后，有关向量存储在工作空间中，在此基础上，观察上述等值线绘制指令的运行效果。
【[cs,h]=contour（X,Y,z,15）;clabel（cs,h,'labelspacing',244）】
2.三元函数可视化:
slice指令
◆观察：
绘制三元函数
的可视化图形。
【clf,x=linspace（-2,2,40）;y=x;z=x;
[X,Y,Z]=meshgrid（x,y,z）;w=X.^2+Y.^2+Z.^2;
slice（X,Y,Z,w,[-1,0,1],[-1,0,1],[-1,0,1]）,colorbar】
3.空间曲线及其运动方向的表现：
plot3和quiver指令
【clf,t=0:
0.1:
1.5;
Vx=2*t;Vy=2*t.^2;Vz=6*t.^3-t.^2;
x=t.^2;y=（2/3）*t.^3;z=（6/4）*t.^4-（1/3）*t.^3;%由速度得到曲线
plot3（x,y,z,'r.-'）,holdon%画飞行轨迹
%算数值梯度，也就是重新计算数值速度矢量，这只是为了编程的方便，不是必须的
图2.12飞机的飞行轨迹与方向
Vx=gradient（x）;Vy=gradient（y）;Vz=gradient（z）;
quiver3（x,y,z,Vx,Vy,Vz）,gridon%画速度矢量图
xlabel（'x'）,ylabel（'y'）,zlabel（'z'）】
2.2应用、思考和练习
2.2.1线性p周期函数
zxy2_3_f.m
【functionf=zxy2_3_f（x）
f=sin（x+cos（x））;】
zxy2_3.m
【clear，clf
a=-8;b=12;n=300;xx=linspace（a,b,n）;
h=zxy2_3_f（xx）;
S
（1）=0;
fori=2:
n
S（i）=S（i-1）+quad（'zxy2_3_f',xx（i-1）,xx（i））;
end
subplot（1,2,1）,plot（xx,S,'k-'）,axis（[a,b,-1.5,9]）
subplot（1,2,2）,plot（xx,[h;zeros（1,length（xx））],'k-'）,axis（[a,b,-1.5,1.5]）】
2.2.2平面截割法和曲面交线的绘制
◆用平行截面法讨论由曲面
构成的马鞍面形状。
zxy2_6.m（平行截割法示例，本程序的绘制两曲面交线方法可以套用）
【clf,a=-20;eps0=1;
[x,y]=meshgrid（-10:
0.2:
10）;%生成平面网格
v=[-1010-1010-100100];%设定空间坐标系的范围
colormap（gray）%将当前的颜色设置为灰色
z1=（x.^2-2*y.^2）+eps;%计算马鞍面函数z1=z1（x,y）
z2=a*ones（size（x））;%计算平面z2=z2（x,y）
r0=abs（z1-z2）<=eps0;
%计算一个和z1同维的函数r0,当abs（z1-z2）<=eps时r0=1；当abs（z1-z2）>eps0时，r0=0。
%可用mesh（x,y,r0）语句观察它的图形，体会它的作用，该方法可以套用。
zz=r0.*z2;xx=r0.*x;yy=r0.*y;%计算截割的双曲线及其对应的坐标
subplot（2,2,2）,%在第2图形窗口绘制双曲线
h1=plot3（xx（r0~=0）,yy（r0~=0）,zz（r0~=0）,'+'）;
set（h1,'markersize',2）,holdon,axis（v）,gridon
subplot（2,2,1）,%在第一图形窗口绘制马鞍面和平面
mesh（x,y,z1）;grid,holdon;mesh（x,y,z2）;
h2=plot3（xx（r0~=0）,yy（r0~=0）,zz（r0~=0）,'.'）;%画出二者的交线
set（h2,'markersize',6）,holdon,axis（v）,
fori=1:
5%以下程序和上面是类似的，通过循环绘制一系列的平面去截割马鞍面
a=70-i*30;%在这里改变截割平面
z2=a*ones（size（x））;r0=abs（z1-z2）<=1;zz=r0.*z2;yy=r0.*y;xx=r0.*x;
subplot（2,2,3）,
mesh（x,y,z1）;grid,holdon;mesh（x,y,z2）;hiddenoff
h2=plot3（xx（r0~=0）,yy（r0~=0）,zz（r0~=0）,'.'）;axis（v）,grid
subplot（2,2,4）,
h4=plot3（xx（r0~=0）,yy（r0~=0）,zz（r0~=0）,'o'）;
set（h4,'markersize',2）,holdon,axis（v）,gridon
end】
2.2.3微分方程的斜率场
◆绘制微分方程
的斜率场，并将解曲线画在图中，观察斜率场和解曲线的关系。
zxy2_5.m（绘制一阶微分方程的斜率场和解曲线）
【clf,clear%清除当前所有图形窗口的图像，清除当前工作空间的内存变量。
a=0;b=4;c=0;d=4;n=15;
[X,Y]=meshgrid（linspace（a,b,n）,linspace（c,d,n））;%生成区域中的网格。
z=X.*Y;%计算斜率函数。
Fx=cos（atan（X.*Y））;Fy=sqrt（1-Fx.^2）;%计算切线斜率矢量。
quiver（X,Y,Fx,Fy,0.5）,holdon,axis（[a,b,c,d]）
%在每一网格点画出相应的斜率矢量，0.5是控制矢量大小的控制参数，可以调整。
[x,y]=ode45（'zxy2_5f',[0,4],0.4）;%求解微分方程。
%zxy2_5f.m是方程相应函数f（x,y）的程序，单独编制；[x0,xs]=[0,4]为求解区间；
%y0=0.4为初始值；输出变量x,y分别为解轨线自变量和因变量数组。
plot（x,y,'r.-'）%画解轨线】
zxy2_5f.m（微分方程的函数子程序）
【functiondy=zxy2_5f（x,y）
dy=x.*y;】
2.2.4颜色控制和渲染及特殊绘图指令
1.地球表面的气温分布（sphere指令）
◆
【[a,b,c]=sphere（40）;t=max（max（abs（c）））-abs（c）;surf（a,b,c,t）;
axis（'equal'）,colormap（'hot'）,shadingflat,colorbar】
2.旋转曲面的生成：
柱面指令cylinder和光照控制指令surfl
◆
【x=0:
0.1:
10;z=x;y=1./（x.^3-2.*x+4）;
[u,v,w]=cylinder（y）;surfl（u,v,w,[45,45]）;
shadinginterp】
3.若干特殊图形
◆运行下面程序，了解各指令的用法和效果。
【x=[1:
10];y=[56348110356];
subplot（2,3,1）,bar（
升级会员 