数学实验与 Matlab 《数学实验与 Matlab》程序 周晓阳 华中科技大学数学系 我将程序按书中的顺序排列,这样便于读者利用。 本书程序均通过了调式。可直接拷贝到命令窗口运行或编制 M文件运行。 如出现问题,可能是中英文标点的缘故(排版有可能使用了 中文标点),请将中文标点换为英文标点试试
数学实验与 Matlab 1 1 《数学实验与 Matlab》程序 周晓阳 华中科技大学数学系 我将程序按书中的顺序排列，这样便于读者利用。 本书程序均通过了调式。可直接拷贝到命令窗口运行或编制 M 文件运行。 如出现问题，可能是中英文标点的缘故（排版有可能使用了 中文标点），请将中文标点换为英文标点试试
数学实验与 Matlab 实验1。矩阵运算与 Matlab命令 11知识要点与背景:日常矩阵及其运算 A=423;132;133;322], %表1-1、表1-2的数据分别写成 矩阵形式 B=[35206045;10155040,20124520 【 C=A*B %矩阵乘法,求各订单所对应的原材料和劳动 【whos %查看 Matlab工作空间中变量及其规模】 12实验与观察:矩阵和 Matlab语言 121向量的生成和运算 x= linspace(0,4*pl,100);%将[0,4x]区间100等分,产生了一个100 维向量 y=sin(x) %计算函数值,产生了一个与x同维的100维函 数向量 y I=sin(x).2 %计算函数向量,注意元素群运算 %以x为横坐标,y为纵坐标画函数的图用不同的线型将函数曲线绘制在一个图 plot(x,y 1.向量的创建 ◆直接输入向量。 【xl=[124],x2=[1,2,1],x3=x1"】 ◆冒号创建向量
数学实验与 Matlab 2 2 实验１．矩阵运算与 Matlab 命令 1.1 知识要点与背景：日常矩阵及其运算 【 A=[4 2 3;1 3 2;1 3 3;3 2 2], % 表1-1、 表1-2的 数据分别写成 矩阵形式 B=[35 20 60 45;10 15 50 40;20 12 45 20] 】 【 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=[1 2 4],x2=[1,2,1],x3=x1' 】 ◆ 冒号创建向量
数学实验与 Matlab x1=3.4:6.7 3.4:2: X3=2.6:-0.8:0】 ◆生成线性等分向量。 x= Flinspace(0,1,5)】 2.向量的运算 r y=sin(x) 【yl=sin(x).^2;y2=exp(-x)*sin(x),】 1.22矩阵创建和运算 1创建矩阵 (1)数值矩阵的创建 ◆直接输入法创建简单矩阵。 【A=[1234;5678;9101112]】 【B=[-1.3,sqrt(3);(1+2)*4/5,sin(5)exp(2),6]】 (2)符号矩阵的创建 sy ms all al2 al3 al4 a21 a22 a23 a24 a31 a32 a33 a34. bl1 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 Al==[all a12 a13 a14: a21 a22 a23 a24; a31 a32 a33 a34 Bl=[bll b12 b13 b14; b21 b22 b23 b24: b31 b32 b33 b34] 2矩阵的运算 【C=A1+B1,D=A1-B1】
数学实验与 Matlab 3 3 【 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=[1 2 3 4; 5 6 7 8; 9 10 11 12] 】 【 B=[-1.3,sqrt(3);(1+2)* 4/5,sin(5);exp(2),6] 】 （2）符号矩阵的创建 ◆ 【 syms a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 … b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 A1=[a11 a12 a13 a14 ;a21 a22 a23 a24; a31 a32 a33 a34], B1=[b11 b12 b13 b14 ;b21 b22 b23 b24; b31 b32 b33 b34] 】 2.矩阵的运算 【 C=A1+B1,D=A1-B1 】
数学实验与 Matlab syr 【C=A1*B1】 ?? Error using == sym/mtimes, Inner matrix dimensions m agree A2=A1(,1:3),B1】 【g11=A2(1,)*B1(:,1)】 H=[123;210;123J,K=[123;210;231 k det=det(K), H inv=inv(H), K inv=K - A=[301;110;014] B=inv(A-2*eye(3)*A, B=(A-2eye(3))\A I 3分块矩阵:矩阵的裁剪、分割、修改与抽取 A=[10112:01-123:30510,23121],wr=[1,3],vc=[1,3] Al=A(vr, vc) %取出A的1、3行和1、3列的交叉处元素构成新矩阵A1】 ◆将上面的矩阵A分为四块,并把它们赋偆铚矩阵B中,观察运行后的结果 【A1l=A(1:2,1:2),A12=A(1:2,3:5),A21=A(3:4,1;2),A22=A(3:4,3:5) B=[All A12; A21 A22]
数学实验与 Matlab 4 4 【 syms c cA=c*A1 】 【 C=A1*B1 】 { ??? Error using ==> sym/mtimes, Inner matrix dimensions must agree. } 【 A2=A1(:,1:3), B1 】 【 G=A2*B1 】 【 g11=A2(1,:)*B1(:,1) 】 【 A, A_trans=A' 】 【 H=[1 2 3 ; 2 1 0 ; 1 2 3 ], K=[1 2 3 ; 2 1 0 ; 2 3 1] h_det=det(H), k_det=det(K),H_inv=inv(H),K_inv=K^-1 】 【 A=[3 0 1; 1 1 0;0 1 4]; B=inv(A-2* eye(3))*A, B=(A-2* eye(3))\A 】 3.分块矩阵：矩阵的裁剪、分割、修改与抽取 （1） 【 A=[1 0 1 1 2;0 1 -1 2 3;3 0 5 1 0;2 3 1 2 1], 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=[A11 A12;A21 A22] 】
数学实验与 Matlab (2)矩阵的修改和提取 【A=[10112;01-123;30510;23121 A(1,)=[00000 ◆观察 B(:[2,4])=[] %删除矩阵B的第2、4列】 (3)矩阵元素的抽取 4生成特殊矩阵 r y1=rand (1, 5), y2=rand (1, 5), rand('seed, 3), xl=rand(1, 5), rand ('seed, 3), x2=rand(1, 5) y 5.常用矩阵函数 6.数据的简单分析 【 rand('seed,D); A Asort=sort(A), Amax=max(A), Asum=sum(A)I 1.23 Matlab工作环境和编程 2Maab的基本设计 13应用、思考与练习
数学实验与 Matlab 5 5 （2）矩阵的修改和提取 ◆ 【 A=[1 0 1 1 2;0 1 -1 2 3;3 0 5 1 0;2 3 1 2 1] A(1,:)=[0 0 0 0 0]; 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.3 Matlab 工作环境和编程 2.Matlab 的基本设计 1.3 应用、思考与练习
数学实验与 Matlab 13.1关系矩阵 132投入产出 133循环比赛的名次 【A=[0110,0011;,0001;,1000], 画矩阵结构图的 gplot指令。 clf,A=[0110:0011:0001;1000];xy=[01:00;-1-0 graphy_ plot(A, xy, 1,0.5),% gplot(A, xy) I 134参考程序 graphy_ plot t function y=graphy plot(A, xy, l,p) %画矩阵的有向结构图。A为邻接矩阵,xy为顶点坐标,1控制参数,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 gplot(A, xy, axis(la bc d, hold on elseif l-=0 U=[;V=[;X=[;Y= length(k)
数学实验与 Matlab 6 6 1.3.1 关系矩阵 1.3.2 投入产出 1.3.3 循环比赛的名次 【 A=[0 1 1 0; 0 0 1 1; 0 0 0 1; 1 0 0 0], e=ones(4,1); c=A* e; s=c' 】 ★ 画矩阵结构图的gplot指令。 ◆（3） 【 clf, A=[0 1 1 0;0 0 1 1;0 0 0 1;1 0 0 0]; xy=[0 1;0 0;-1 –0.5;1 –0.5]; graphy_plot(A,xy,1,0.5), % gplot(A,xy) 】 1.3.4 参考程序 graphy_plot.m 【 function y=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; if l=0 gplot(A,xy),axis([a b c d]),hold on, elseif l~=0 U=[];V=[];X=[];Y=[]; n=length(A(:,1)) ; for i=1:n k=find(A(i,:)~=0); m=length(k);
数学实验与 Matlab if(m~=0) u(1)=(xy(k(),1)-xy(i,1))v(1)=(xy(k(j),2)-xy(i,2)) u(2)=ep X=[[xy(i,1)xy(k(j),1)]X],Y=[xy(i,2)xy(k(j),2)]:Y text(xy(i, 1),xy (i, 2),'\bullet\leftarrowlfontsize( 16)itV) m2str(i)]): hold gplot(A, xy), axis([a bc d)), hold on h=quiver(x,Y, U, V, p); set(h, color,red); hold o plot(xy (: 1), xy(: 2), k !,markersize, 12), hold end hold off 实验2.函数的可视化与 Matlab作 21实验与观察:函数的可视化 211 Matlab二维绘图命令 1周期函数与线性p一周期函数 观察 r clf, x-linspace(0, 8*pi, 100) F=inline(sin(x+cos(x+sin(x)))) yI=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))
数学实验与 Matlab 7 7 if(m~=0) for j=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)]); hold on, end end gplot(A,xy),axis([a b c d]),hold on, h=quiver(X,Y,U,V,p);set(h,'color','red');hold on, plot(xy(:,1),xy(:,2),'k.','markersize', 12),hold on, end , hold off 】 实验２．函数的可视化与 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) 】
数学实验与 Matlab 2plot指令:绘制直角坐标的二维曲线 3.图形的属性设置和屏幕控制 h=plot([0:.0.1:2*pil,sin([0:0.1:2*pi]);grid set(h, ' Line width, 5, color,'red) set(gca, 'Grid LineStyle, '-' 'fontsize, 16)I ◆设置y坐标的度并加以说明,并改变字体的大小 h=plot(0:0.1:2*pi],sin([0:0.1:2*pi])grid set(gca,'ytick,[-1-05005 1), set(gca, 'yticklabel,'alblcldJe), 4.文字标注指令 plot(x,yl, ,b, x,y2, k-) se(gca," fontsize',l5, ontname', times new roman),%设置轴对象的 体为tme % New roman,字体的大小为15 title( lit( Peroid and linear peroid function)); %加标题,注意文字用单引号 加上 %斜杠'后可输入不同的设置,例如it{…}表示花括号里的文字为斜体;如果有多项设 %则可用\...连续输入 xlabel('x from 0 to 8*pi it(t)\); ylabel('\it(y),) %说明坐标轴 text(x(49),y 1(50)-0.4, \fontsize( 15)\bullet\leftarrow The period function itf(x)}"); %在坐标(x(49)y1(50)-0.4)处作文字说明,各项设置用"\"隔开。 % fontsize{l5} bullet\leftarrow的意义依次是:1字体大小=15\画圆点 左箭头 text(x(14),y2(50)+1, \fontsize( 15) The linear function fitg(x) rightarrow \bullet) %与上一语句类似,用右箭头
数学实验与 Matlab 8 8 2. plot 指令：绘制直角坐标的二维曲线 3. 图形的属性设置和屏幕控制 【 h=plot([0:0.1:2* pi],sin([0:0.1:2* pi])); grid on 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]));grid on， set(gca,'ytick',[-1 -0.5 0 0.5 1]), 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','times New Roman'), %设置轴对象的 字体为 times % New Roman，字体的大小为 15 title(' \it{Peroid and linear peroid function} '); %加标题，注意文字用单引号 ' ' 加上 %斜杠'\'后可输入不同的设置，例如it{…}表示花括号里的文字为斜体；如果有多项设 置， %则可用\…\…\…连续输入。 xlabel('x from 0 to 8* pi it{t}\'); ylabel('\it{y}'); %说明坐标轴 text(x(49),y1(50)-0.4,'\fontsize{15}\bullet\leftarrowThe period function {\itf(x)}'); %在坐标(x(49),y1(50)-0.4）处作文字说明， 各项设置用"\"隔开。 %\fontsize{15}\bullet\leftarrow的意义依次是：\字体大小=15 \ 画圆点 \ 左箭头 text(x(14),y2(50)+1,'\fontsize{15}The linear period function {\itg(x)}\rightarrow \bullet') %与上一语句类似，用右箭头 】
数学实验与 Matlab ◆观察指令 legend和mum2s的用 Peroid and linear peroid function 法:在同张图上画出y=sn()和y=3e0, The linear period functiong(x 这里t∈[O,3],并进行适当的标主 zxy2_2.m If t=0:0.1:3*pi; alpha=0:0.1:3*pi 101520253 plot(t, sin(t), 'r-); hold or 图2.5文字标注 plot(alpha, 3*exp(-05*alpha), ' k: ' set(gca, 'fontsize, 15, ' fontname, ' times New Roman) xlabel (\it(t(deg),'); ylabel('\it(magnitude); title('lit(sine wave and f\it(Ae))t-alphaflitt wave); %i &r \alphal 意义 text(6, sin(6), \fontsize(15) The Value \(sin(t))at fitt =6\rightarrowlbullet', 'HorizontalAlignment', 'right,) %上面的语句是一整行,如果要写成两行,必须使用续行号…,例如要在 bullet %后换行,需写“ bullet',…”后才能换行。 %' Horizontalalignment';righ!'表示箭头所指的曲线对象在文字的右边。 text(2, 3*exp(-0 5*2), '\fontsize 15)\bullet\leftarrow The Value of ); %num2str的用法:[' string',num2str,' string2],注意方括号的使用 % legend(" Itsin(t)';{ \itAe}^{- \alphat}")%请结合图形观察此命令的使 运行结果如图26所示 5.图形窗口的创建和分割 观察 clf, b=2+pi; x=linspace(0, b, 50)
数学实验与 Matlab 9 9 ◆ 观察指令 legend 和 num2str 的用 法：在同一张图上画出 t y t y e 0.5 sin( ) 3 − = 和 = ， 这里 t [0,3], 并进行适当的标注。 zxy2_2.m 【 clf, t=0:0.1:3* pi;alpha=0:0.1:3* pi; plot(t,sin(t),'r-');hold on; plot(alpha,3* exp(-0.5* alpha),'k:'); set(gca,'fontsize',15,'fontname','times New Roman'), xlabel('\it{t(deg)}');ylabel('\it{magnitude}'); title(' \it{sine wave and { \it{Ae}}^{-\alpha{\itt}}wave}'); %注 意\alpha的 意 义 text(6,sin(6),'\fontsize{15}The Value \it{sin(t)} at {\itt}=6\rightarrow\bullet', 'HorizontalAlignment','right'), %上面的语句是一整行，如果要写成两行，必须使用续行号 … ，例如要在 “ bullet',” %后换行，需写“bullet', …”后才能换行。 % 'HorizontalAlignment','right' 表示箭头所指的曲线对象在 文字的右边。 text(2,3* exp(-0.5* 2),['\fontsize{15}\bullet\leftarrow The Value of \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); 0 5 10 15 20 25 30 -1 0 1 2 3 4 5 6 Peroid and linear peroid function x y ?? The period function f(x) The linear period function g(x)? ? 图 2.5 文字标注
数学实验与 Matlab for k=1: 9 subplot(3, 3, k), plot(x, y), axis([0, 2*pi, -1, ID) 212多元函数的可视化与空间解析几何(三维图形) 本节通过高等数学的几个例子观察 Matlab的三维绘图功能和技巧 1.绘制二元函数 ◆观察:绘制κ=∫(x,y)=(1-x)h(x-y)的图象,作定义域的栽剪 ◆(1)观察 meshgrid指令的效果。 【a=-0.98b=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、 mesha、surf ◆(2)做團数的定义域裁剪,瘀察上述三维绘图指令的效果。 程序zxy24m clear clf a=-1;b=1;c=-15;d=15;n=20;epsl=0.01 x=linspace(a, b, n); y=linspace(c, d, n) [X,Y=meshgrid(x,y) %计算函数值z,并作 定义域裁剪 if(1-X(1j))<epsl|x(i)-Y(1,j)<epsl%i语句这样用 (1,j)=NaN; %作定义域 裁剪,定义域以外的函数值为NaN
数学实验与 Matlab 10 10 for k =1:9 y=sin(k* x); subplot(3,3,k),plot(x,y),axis([0,2* pi,-1,1]) end 】 2.1.2 多元函数的可视化与空间解析几何（三维图形） 本节通过高等数学的几个例子观察 Matlab 的三维绘图功能和技巧。 1. 绘制二元函数 ◆ 观察：绘制 ( , ) (1 ) ln( ) 2 1 z = f x y = − x x − y − 的图象，作定义域的裁剪。 ◆ （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); for i=1:n %计算函数值z ，并 作 定义域裁剪 for j=1:n if (1-X(i,j))<eps1|X(i,j)-Y(i,j)<eps1 %if语句这样用 z(i,j)=NaN; %作定义域 裁剪，定义域以外的函数值为NaN