数学实验与matlab 程序解析资料全编.docx
《数学实验与matlab 程序解析资料全编.docx》由会员分享,可在线阅读,更多相关《数学实验与matlab 程序解析资料全编.docx(85页珍藏版)》请在冰豆网上搜索。
![数学实验与matlab 程序解析资料全编.docx](https://file1.bdocx.com/fileroot1/2023-5/19/a822c174-7821-4412-a659-e6e8178623af/a822c174-7821-4412-a659-e6e8178623af1.gif)
数学实验与matlab程序解析资料全编
《数学实验与Matlab》程序
将程序按书中的顺序排列,本书程序均通过了调式。
可直接拷贝到命令窗口运行或编制M文件运行。
如出现问题,可能是中英文标点的缘故(排版有可能使用了中文标点),请将中文标点换为英文标点试试。
实验1.矩阵运算与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】
实验2.函数的可视化与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(