-matlab在理工课程中的应用

1、MATLAB在理工课程中的应用第五章 在高等数学中的应用5.1 函数、极限和导数一、单变量函数的计算和绘图例5-1-1设 要求以0.01秒为间隔，求出y的151个点，并求出其导数的值和曲线。解： 建模 直接用程序文件编程的方法 编成函数文件，由主程序调用的方法t=0:.01:1.5; %设定自变量数组tw=4*sqrt(3); %固定频率y=sqrt(3)/2*exp(-4*t).*sin(w*t+pi/3);%注意用数组运算式subplot(2,1,1),plot(t,y),grid %绘制曲线并加上坐标网格title(绘图示例),xlabel(时间t),ylabel(y(t) % 加标注%

2、 求导数并绘制导数曲线，注意数组求导数后其长度减少1Dy=diff(y);subplot(2,1,2),plot(t(1:length(t)-1),Dy),gridylabel(Dy(t) % 加标注编成函数文件：%主程序(ex511b)dt=0.01;t=0:dt:1.5;w=4*sqrt(3);y=ex511bf(t,w);Dy=diff(y)/dt;subplot(2,1,1),plot(t,y),grid %绘制曲线并加上坐标网格subplot(2,1,2),plot(t(1:length(t)-1),Dy),gridylabel(Dy(t) % 加标注%函数文件ex511bffunc

3、tion xvalues=ex511bf(tvalues,w)xvalues=sqrt(3)/2*exp(-4*tvalues).*sin(w*tvalues+pi/3); 从本例来看，第二种方法似乎更麻烦一些，但它具有模块化的特点。当程序中要反复多次调用此函数，而且输入不同的自变量时，利用函数文件可大大简化编程。例5-1-2绘制极坐标系下曲线并讨论参数a,b,n的影响。为了便于比较，编制一个能分别画两个图形的程序，采用for循环。 可从中看出用循环的技巧%绘制极坐标系下曲线theta=0:0.1:2*pi; %产生极角向量for i=1:2 a(i)=input(a=);b(i)=input

4、(b=);n(i)=input(n=) rho(i,:)=a(i)*cos(b(i)+n(i)*theta); %极坐标方程 subplot(1,2,i),polar(theta,rho(i,:);%极坐标系绘图enda=2 b=pi/4 n=2(4叶玫瑰线)a=2 b=0 n=3(3叶玫瑰线)二、参变方程表示的函数的计算和绘图例5-1-3摆线的绘制。当圆轮在平面上滚动时，其圆周上任一点所画的轨迹称为摆线。 这一点在圆内 内摆线 这一点在圆外 外摆线rtARr=1XYOr为轮半径，R为点半径t=0:0.1:10;r=input(r=),R=input(R=)x=r*t-R*sin(t);y=r

5、-R*cos(t);plot(x,y),hold on,axis equal %x,y坐标保持等比例显示摆线的绘制例5-1-4三次抛物线的方程为 试探讨参数a和c对其图形的影响。三、曲线族的绘制x=-2:0.1:2; %给定x数组，确定范围及取点密度subplot(1,2,1) %分两个画面绘图for c=-3:3 plot(x,x.3+c*x),hold on,end,grid %a=1 取不同的caxis(equal),axis(-2 2 -3 3) %x,y坐标等比例并确定范围subplot(1,2,2)for a=-3:3 plot(x,a*x.3+x),hold on,end,gri

6、d %c=1 取不同的aaxis equal ,axis(-2 2 -3 3)c和a取不同值时的曲线族5.2 空间解析几何一、曲面方程例5-2-1二次曲面的方程如下 要求讨论参数a,b,c对其形状的影响，并画出其图形。a=input(a=);b=input(b=);c=input(c=);d=input(d=);N=input(N=); %输入参数，N为网格线数目xgrid=linspace(-abs(a),abs(a),N); %建立x网格坐标ygrid=linspace(-abs(b),abs(b),N); %建立y网格坐标x,y=meshgrid(xgrid,ygrid) %确定NXN个

7、点的x,y网格坐标z=c*sqrt(d-y.*y/b/b-x.*x/a/a);u=1; %u=1,表示z要取正负值z1=real(z); %取z的实部z1for k=2:N-1 %以下7行程序的作用是取消z中含虚数的点 for j=2:N-1 if imag(z(k,j)=0 z1(k,j)=0;end if all(imag(z(k-1:k+1,j-1:j+1)=0 z1(k,j)=NaN;end endendsurf(x,y,z1),holdif u=1 z2=-z1;surf(x,y,z2) %u=1时加画负半面，并加坐标轴 axis(-abs(a),abs(a),-abs(b),abs

8、(b),-abs(c),abs(c);endxlabel(x),ylabel(y),zlabel(z)hold offa=5 b=4 c=3 d=1 N=20a=5i b=4 c=3 d=1 N=20a=5i b=4i c=3 d=1 N=20二、空间两曲线的交线例5-2-2列出求空间任意曲面的交线的程序 两空间曲面联立起来就形成一个空间的曲线的方程 mesh语句 要显示其交点，须找到各个交点。由于数值计算是离散点，难以找到完全重合的点 采用设置门限的方法，两曲线差值小于设定门限，就认为是交点x,y=meshgrid(-2:.1:2); %确定计算和绘图的定义域网格z1=x.*x-2*y.*y

9、; %第一个曲面方程z2=2*x-3*y; %第二个曲面方程(平面)mesh(x,y,z1);hold;mesh(x,y,z2); %在一个图上同时画两个图r0=(abs(z1-z2)=.1; %求二个曲面z坐标差小于0.1的网格矩阵zz=r0.*z1;yy=r0.*y;xx=r0.*x; %求这些网格上的坐标值，即交线坐标plot3(xx(r0=0),yy(r0=0),zz(r0=0),*);%画出这些点colormap(gray),hold off两曲面的交线 如果想改变曲面方程，可以改动第二行和第三行，但这样的程序不是通用的。最好是程序运行时能向用户提问，允许用户输入曲面方程 用到字符串

10、功能和eval命令 S1=input(输入第一个方程,s); 在原来的z1语句改为 z1=eval(s1); 给出定义域和间隔 x,y=meshgrid(xmin:dx:xmax，ymin:dy:ymax);xmin=input(xmin=);dx=input(dx=);xmax=input(xmax=)ymin=input(ymin=);dy=input(dy=);ymax=input(ymax=)x,y=meshgrid(xmin:dx:xmax,ymin:dy:ymax); %确定计算和绘图的定义域网格s1=input(输入第一个方程:,s); %输入第一个曲面方程z1=eval(s1)

11、;s2=input(输入第二个方程:,s); %输入第二个曲面方程z2=eval(s2);mesh(x,y,z1);hold;mesh(x,y,z2); %在一个图上同时画两个图r0=(abs(z1-z2)=.1; %求二个曲面z坐标差小于0.1的网格矩阵zz=r0.*z1;yy=r0.*y;xx=r0.*x; %求这些网格上的坐标值，即交线坐标plot3(xx(r0=0),yy(r0=0),zz(r0=0),*);%画出这些点colormap(gray),hold off5.3 数列和极限5.4 数值方法和数值积分一、任意非线性方程f(x)=0的解例5-4-1用切线法求方程的近似数值解即求任

12、意曲线y=f(x)过零点的问题fzero先用fplot函数得出曲线，大致了解一下曲线形状fplot(x.3+10*x.2-2*sin(x)-50,-12,5)%子程序function y=ex541f(x)y=x.3+10*x.2-2*sin(x)-50;%主程序x0=input(x0=);fplot(x.3+10*x.2-2*sin(x)-50,-12,5),hold on,grid onx=fzero(ex541f,x0)plot(x, 0,*r);hold on;x0=-10 x = -9.4384x0=3x = 2.0707x0=-3x = -2.564二、数值定积分trapz函数三、

13、多重积分例5-4-3计算二重积分 积分区域为由x=1,y=x及y=0所围成的闭合区域xyy=x01 建模：积分区域如图所示clear,format compactfill(0,1,1,0,0,0,1,0,y),hold %画出积分区域fill(0.55,0.6,0.6,0.55,0.55,0,0,0.6,0.55,0,r)%画出单元条dx=input(步长dx=);dy=dx;x=0:dx:1;lx=length(x);for k=1:lx x1=(k-1)*dx; y1=0:dy:x1; f=x1.2+y1.2; s1(k)=trapz(f)*dy;ends=trapz(s1)*dx步长dx

14、=0.1s = 0.3375步长dx=0.01s = 0.3334例5-4-4计算三重积分 积分区域为由x=1,y=x，z=xy及z=0所围成的闭合区域 建模： 先画出积分区域三重积分区域图dx=input(步长dx=);dy=dx;dz=dx;x=0:dx:1;lx=length(x);for k=1:lx x1=(k-1)*dx; y1=0:dy:x1; for j=1:length(y) y1=(j-1)*dy; z1=0:dz:x1*y1; f=x1.*y1.2.*z1.3; s1(j)=trapz(f)*dz; ends2(k)=trapz(s1)*dyends=trapz(s2)*

15、dx步长dx=0.1 s = 4.2081步长dx=0.01 s = 0.0333四、微分方程和数值积分例5-4-5用数值积分法求解下列微分方程设初始时间t0=0；终止时间tf=3;初始条件y(0)=0,dy(t)/dt=0。 建模： 先将方程化成两个一阶微分方程组，设写成矩阵形式为clf,t0=0;tf=3*pi;x0t=0;0;t,x=ode23(ex545f,t0,tf,x0t)y=x(:,1); %数值积分解for I=1:length(t); y2(I)=(1+2/(pi2)*(1-cos(t(I)-t(I)2/(pi2);%解析解endu=1-(t.2)/(pi2);clf,plo

16、t(t,y,-,t,u,+,t,y2,o)legend(数值积分解,输入量,解析解)function xdot=ex545f(t,x)u=1-(t.2)/(pi2);xdot=0 1;-1 0*x+0 1*u;主程序：函数程序：5.5 线性代数矩阵的行列式det矩阵的秩rank矩阵的迹trace第六章 在普通物理中的应用6.1 物理数据处理 例6-1-2 写出一个程序，能把用户输入的长度单位在厘米，米，千米，英寸，英尺，英里，市尺，市里之间任意转换建模： 第一步把输入量变换为米 第二步把米变换为输出单位 把变换常数直接表示成一个数组，选择单位的序号也就成了数组的下标clear all; dis

17、p( 长度单位换算程序) fprintf(长度单位: n); %选择输入输出的单位fprintf( 1) 厘米 2) 米 3) 千米 4) 英寸 n);fprintf( 5) 英尺 6) 英哩 7) 市尺 8) 市里 n);InUnits = input(选择输入单位编号: );OutUnits = input(选择输出单位编号: );% 令各种单位对米的变换常数数组 ToMeter ToMeter = 0.01, 1.00, 1000.0, 0.0254, 0.3048, 1609.3, 1/3, 500 ;FrmMeter= 1./ ToMeter;%反变换常数数组FrmMeter为ToM

18、eter数组的倒数Value = input(输入待变换的值(0为退出): ), while( Value = 0 ) ValueinM = Value*ToMeter(InUnits); %把输入值变为米 NewValue = ValueinM*FrmMeter(OutUnits); % 把米变为输出单位 fprintf(变换后的值是 %g n,NewValue); % 打印变换后的值 Value = input(输入待变换的值(0为退出): ); % 提问下个输入值end 6.2 力学基础例6-2-1 设目标相对于设点的高度为yf，给定初速，试计算物体在真空中飞行的时间和距离建模：求落点时

19、间tf时，需要求解一个二次线性方程由解出t，它就是落地时间tf,tf有两个解，我们得到一个有效解（大值），再求clear; y0 = 0; x0 = 0; % 初始位置vMag = input(输入初始速度 (m/s):(书上为50) ); % 输入初始速度vDir = input( 输入初速方向(度):(书上为40或50) );yf = input(输入目标高度（米）:(书上为8) ); % 输入目标高度vx0 = vMag*cos(vDir* (pi/180); % 计算x,y方向的初始速度vy0 = vMag*sin(vDir* (pi/180); % wy = -9.81; wx =

20、0; % 重力加速度 (m/s2)tf=roots(wy/2,vy0,y0-yf); % 解方程wy*t2/2+vy0*t+y0=yh，计算落点tftf=max(tf); % 去除落点时间tf中的庸解t=0:0.1:tf,tf; % 设定时间数组，因tf不大可能被0.1整除，必须加一个tf点y = y0 + vy0*t + wy*t.2/2; % 计算轨迹x = x0 + vx0*t + wx*t.2/2;xf = max(x),plot(x,y),grid % 计算射程，画出轨迹set(gcf,color,w) % 置图形背景色为白色 例621 无阻力抛物体的飞行轨迹6.4 静电场 例6-4

21、-1计算平面上N个电荷之间的库仑引力。 库仑定律是描述真空中两个静止点电荷之间相互作用的实验定律clear all;N=input(输入电荷数目N=:);for ic=1:N fprintf(-n 对电荷 # %g/n,ic); rc=input(输入电荷位置x y(米):); x(ic)=rc(1); % 电荷ic的x坐标 y(ic)=rc(2); q(ic)=input(输入电荷量(库仑):);endE0=8.85e-12;C0=1/(4*pi*E0)for ic=1:N Fx=0.0;Fy=0.0; %初始化力 for jc=1:N if(ic=jc) xij=x(ic)-x(jc);y

22、ij=y(ic)-y(jc); Rij=sqrt(xij2+yij2); Fx=Fx+C0*q(ic)*q(jc)*xij/Rij3; Fy=Fy+C0*q(ic)*q(jc)*yij/Rij3; endendfprintf(其它电荷作用在电荷 # %g上的合力为:n,ic);fprintf(x-分量:%g Nn,Fx);fprintf(y-分量:%g Nn,Fy);end结果：输入电荷数目N=:2- 对电荷 # 1/n输入电荷位置x y(米):1 1输入电荷量(库仑):0.5- 对电荷 # 2/n输入电荷位置x y(米):2 3输入电荷量(库仑):0.7C0 = 8.9918e+009其它电

23、荷作用在电荷 # 1上的合力为:x-分量:-2.81488e+008 Ny-分量:-5.62976e+008 N其它电荷作用在电荷 # 2上的合力为:x-分量:2.81488e+008 Ny-分量:5.62976e+008 N例6-4-3 由电位的表示式计算电场并画出等电位线和电场方向建模 如果已知空间的电位分布 则空间的电场等于电位场的负梯度其中 分别为x,y,z三个方向的单位向量MATLAB设有gradient函数，它是靠数值微分，因此空间观测点应该取得密一些，已获得较高的精度fprintf(输入电位分布方程 V(x,y) n);fprintf(例如: log(x.2 + y.2) n);

24、V = input(: ,s); % 读入字符串 V(x,y)NGrid = 20; % 绘图的网格线数xMax = 5; % 绘图区从 x= -xMax 到 x= xMax yMax=5; xPlot = linspace( -xMax, xMax, NGrid); % 绘图取的x值x,y=meshgrid(xPlot);VPlot=eval(V);ExPlot, EyPlot = gradient(-VPlot); % 电场是电位的负梯度clf; subplot(1,2,1),meshc(VPlot); % 画含等高线的三维曲面set(gcf,color,w) % 置图形背景色为白色xla

25、bel(x); ylabel(y); zlabel(电位);% 规定等高线图的范围及比例subplot(1,2,2), axis(-xMax xMax -yMax yMax); % 建立第二子图cs = contour(x,y,VPlot); % 画等高线clabel(cs); hold on; % 在等高线图上加上编号等高线值% 在等高线图上加上电场方向quiver(x,y,ExPlot,EyPlot); % 画电场 E 的箭头图xlabel(x); ylabel(y);hold off; 电位三维立体图 等位线及电场分布图6.5 稳恒磁场例6-5-1用毕奥萨伐定律计算电流环产生的磁场zId

26、ldBdBPrI0真空的磁导率0410-7H/m。任一电流元Idl在空间任一点P处所产生的磁感应强度为clear all; % 清工作空间及变量初始化mu0 = 4*pi*1e-7; % 真空导磁率 (T*m/A)I0 = 5.0; % 环中电流(A)Rh = input(输入环半径Rh(m):(书上取2) );C0 = mu0/(4*pi) * I0; % 归并常数xMax = 3; yMax = 3; % 规定图的范围NGx = 21; NGy = 21; % 规定观测点网格线数x=linspace(-xMax, xMax, NGx); % 确定观测点的x,y坐标数组y=linspace(

27、-yMax, yMax, NGy);Nh = 20; % 电流环分段数 % 计算每段的端点,环在x=0平面上,其坐标x1,x2均为零theta0 = linspace(0,2*pi, Nh+1); % 环的圆周角分段theta1 = theta0(1:Nh);y1 = Rh*cos(theta1); % 环各段向量的起点坐标y1,z1z1 = Rh*sin(theta1);theta2 = theta0(2:Nh+1);y2 = Rh*cos(theta2); % 环各段向量的终点坐标y2,z2z2 = Rh*sin(theta2); dlx = 0; % 计算环各段向量dl的三个分量dly

28、= y2-y1;dlz = z2-z1; xc = 0; % 计算环各段向量中点的三个坐标分量yc = (y2+y1)/2;zc = (z2+z1)/2; % 循环计算各网格点上的B(x,y) 值for i=1:NGy for j=1:NGx % 对yz平面内的电流环分段作元素群运算,先算环上某段与观测点之间的向量r rx = x(j) - xc; ry = y(i) - yc; rz = -zc; % 观测点在z=0平面上 r3 = sqrt(rx.2 + ry.2 + rz.2).3; % 计算r3 dlXr_x = dly.*rz - dlz.*ry; % 计算叉乘积dl X r的 x

29、和 y 分量 dlXr_y = dlz.*rx - dlx.*rz; Bx(i,j) = sum(C0*dlXr_x./r3); % 把磁场各段的x 和 y 分量累加 By(i,j) = sum(C0*dlXr_y./r3); endendclf; quiver(x,y,Bx,By); % 用quiver 画磁场向量图set(gcf,color,w) % 置图形背景色为白色hold on; plot(0,Rh,bo);plot(0,-Rh,rx); % 在图上画出电流环xlabel(x); ylabel(y);hold off; 电流环产生的磁场分布图6.6 振动与波例6-6-1振动的合成及拍

30、频现象 分别输入两个正弦波的振幅、相位及频率，观察其合成结果，特别是两个信号的频率接近时产生的拍频现象建模： 两个振动相加：两个同方向的振动分别为：用三角函数关系，可求出振动方程t=0:0.001:10;%10秒分10000个点a1=input(振幅1);w1=input(相位1);a2=input(振幅2);w2=input(相位2);y1=a1*sin(w1*t);y2=a2*sin(w2*t);y=y1+y2;subplot(3,1,1),plot(t,y1),ylabel(y1)subplot(3,1,2),plot(t,y2),ylabel(y2)subplot(3,1,3),plo

31、t(t,y),ylabel(y),xlabel(t)pause,sound(y1);pause(2),sound(y2);pause(2),sound(y),pausesubplot(1,1,1)振幅11.2相位1300振幅21.5相位22906.7 光学例6-7-1两点（双缝）光干涉图案 单色光通过两个窄缝射向屏幕，相当于位置不同的两个同频同相光源向屏幕照射的叠和。由于到达屏幕各点的距离不同，引起相位差，有的点加强，有的点抵消，造成干涉现象。yysOzL2L1S2S1光缝屏幕双缝干涉示意图振幅光强两光源到屏幕的距离光程差相位差clcLambda = input(输入光的波长(单位为 nm):

32、 (书上取500) );Lambda = Lambda * 1e-9; % 将nm换 变为 md = input(输入两个缝的间距 (单位为 mm): (书上取2) );d = d * 0.001; % 将mm 变换为 mZ = input(输入缝到屏的距离 (单位为 m): (书上取1) );yMax = 5*Lambda*Z/d; xs = yMax; % 设定图案的y,x 向范围Ny=101;ys = linspace(-yMax,yMax,Ny); % y方向分成101点for i=1:Ny % 对屏上全部点进行循环计算 % 计算第一和第二个光源到屏上各点的距离 L1 = sqrt(y

33、s(i)-d/2).2 + Z2 ); L2 = sqrt(ys(i)+d/2).2 + Z2 ); Phi = 2*pi*(L2-L1)/Lambda; % 从距离差计算相位差 B(i,:) = 4*cos(Phi/2).2; % 计算该点光强（设两束光强相同)end% 在屏上画出图象% clf; figure(gcf); % 清图形窗，将它移到前面NCLevels = 255; % 确定用的灰度等级% 定标：使最大光强(4.0)对应于最大灰度级（白色）Br = (B/4.0) * NCLevels;subplot(1,4,1),image(xs,ys,Br); % 画图象colormap(

34、gray(NCLevels); % 用灰度级颜色图subplot(1,4,2), plot(B(:),ys) % 画出沿y 向的光强变化曲线波长(m)0.0000005光缝距离(m)0.002光栅到屏幕距离(m)1双缝干涉条纹及光强分布第七章 在力学、机械中的应用例7-1-4四连杆机构如图所示：输入杆L1的转角 求输出杆L3转角随时间的变化规律，并求其角速度和角加速度。建模：四连杆机构的运动方程X和Y方向的长度关系确定为（1）（2）从上述两个方程中削去 ，便可划成一个只包括 和 的方程，给定 ，可求出满足此方程的由（2）式得 将（1）中的 代以 得出求得 、 和 后，就可以根据杆1的角速度求出

35、杆3的角速度为了求能使f=0的 值，利用fzero函数，把f=f（ ） 定义为函数文件ex714f.m 利用global命令定义了杆件的长度参数 function y=ex714f(x)global L0 L1 L2 L3 th1y=L1.*cos(th1)+L2*sqrt(1-(L3*sin(x)-L1*sin(th1).2/L2/L2)-L3*cos(x)-L0; 主程序ex714a.mglobal L0 L1 L2 L3 th1 L0=20;L1=8;L2=25; L3=20; % 输入基线及三根杆的长度L1,L2,L3theta1=input(当前角theta1=(书上取0) );w1

36、=input(杆1角速度w1=(书上取100) );theta3=input(对应于theta1的theta3近似值=(书上取1) );th1=theta1;theta3=fzero(ex714f,theta3); % 求初始输出theta3theta2 = asin( L3*sin(theta3)- L1*sin(theta1)/L2); % 计算三角关系w3 = L1*w1*cos(pi/2-theta1+theta2)/ (L3*cos(theta3-pi/2-theta2) ex714b.m 计算输入输出杆件的角位置关系和输出角速度global L0 L1 L2 L3 th1 L0=2

37、0;L1=8;L2=25; L3=20;% 输入基线及三根杆的长度L1,L2,L3w1=input(杆1角速度w1= (书上取100) 1/秒);theta1=linspace(0,2*pi,181); %杆1每圈分为180份，间隔2度。theta3=input(对应于theta1最小值处的theta3近似值= (书上取1) 弧度);dt = 2*pi/180/w1; % 杆1转2度对应的时间增量th1=theta1(1);theta3(1)=fzero(ex714f,theta3); %求初始输出theta3for i=2:181 th1=theta1(i); theta3(i)=fzero

38、(ex714f,theta3(i-1); % 调用ex714f函数逐次求输出theta3endsubplot(1,2,1),plot(theta1,theta3);xlabel(theta1),ylabel(theta3),grid % 画曲线set(gcf,color,w) % 置图形背景色为白色w3 = diff(theta3)/dt; % 求杆3的角速度,注意求导数后数组长度小一subplot(1,2,2),plot(theta1(2:length(theta1),w3);grid 四连杆机构的输入输出角位置关系和输出角速度% 四连杆运动的分析计算:求全程运动并作动画global L0

39、L1 L2 L3 th1 L0=20;L1=8;L2=25; L3=20; % 输入基线及三根杆的长度L1,L2,L3w1=100; % input(杆1角速度w1= (书上取100) 1/秒);theta1=linspace(0,2*pi,181); % 杆1每圈分为180份，间隔2度。theta3=1; % input(对应于theta1最小值处的theta3近似值= (书上取1) 弧度);dt = 2*pi/180/w1; % 杆1转2度对应的时间增量th1=theta1(1);theta3(1)=fzero(ex714f,theta3); % 求初始输出theta3for i=2:18

40、1 th1=theta1(i); theta3(i)=fzero(ex714f,theta3(i-1); % 调用ex714f函数逐次求输出theta3endfigure(1)set(gcf,color,w) % 置图形背景色为白色subplot(1,2,1),plot(theta1,theta3);xlabel(theta1),ylabel(theta3),grid %画曲线w3 = diff(theta3)/dt; % 求杆3的角速度,注意求导数后数组长度小一subplot(1,2,2),plot(theta1(2:length(theta1),w3);gridpause,figure(2

41、),subplot(1,1,1),axis equal % 设定第二张图，使纵横比例相同axis(0,50,-20,20),axis offx0=10;y0=0; % 左支点坐标x1=x0+L1*cos(theta1);y1=L1*sin(theta1); % 杆1右端点坐标x2=x0+L0+L3*cos(theta3);y2=L3*sin(theta3); % 杆3左端点坐标，即杆二右端点坐标x3=x0+L0;y3=0; % 杆3右端点坐标line(5,45,0,0) % 画出基准线h1=line(x0,x1(1),x2(1),x3(1),y0,y1(1),y2(1),y3,Linewidt

42、h,3); % 将四个端点相联set(h1,erasemode,xor); % 设定更新数据的模式，把前一组数据的图形擦除set(gcf,color,w) % 置图形背景色为白色for i1=1:1000 i=mod(i1,180)+1; % 在前面算出的180度范围内的数据循环取值mod(a,b)就 是求的是a除以b的余数 set(h1,Xdata,x0,x1(i),x2(i),x3,Ydata,y0,y1(i),y2(i),y3) % 更新杆端数据 fft(rand(214,1); %延时， 等待一段时间，减慢更新速度 drawnow % 及时显示end7.2材料力学例7-2-2 长为L的

43、悬臂梁，左端固定，在离固定端L1处施加力P，求它的转角和挠度。设梁的弹性模量E=200*109N/m2,惯性矩I=2*10-4m4建模：从转矩求转角要经过一次不定积分，从转角求挠度又要经过一次不定积分。在matlab中可通过cumsum作近似的不定积分，更为精确的函数，还可用cumtrapz函数转角弯矩方程挠度（0 x L1）（L1x L）程序ex722.mclearL=2; P=2000; L1=1.5; % 给出常数E = 200e9; I=2e-5;x = linspace(0,L,101); dx=L/100; % 将x分100段,步长为L/100n1=L1/dx+1; % 确定x=L

44、1处对应的下标M1 = -P*( L1-x(1:n1); % 第一段弯矩赋值M2 = zeros(1,101-n1); % 第二段弯矩赋值（全为零）M = M1,M2; % 全梁的弯矩A = cumsum(M)*dx/(E*I); % 对弯矩积分求转角Y = cumsum(A)*dx; % 对转角积分求挠度subplot(3,1,1),plot(x,M),grid % 绘弯矩图subplot(3,1,2),plot(x,A),grid % 绘转角图subplot(3,1,3),plot(x,Y),grid % 绘挠度图7.3 机械振动例731 分析单自由度阻尼系统的阻尼系数对其固有振动模态的影

45、响建模质量为m、阻尼系数为c、弹性系数为k的单自由度系统自由衰减振动时，其运动微分方程为可改写为式中：系统固有频率； 阻尼比 小阻尼（1）时，微分方程式的解可写为 式中：、由初始条件确定的积分常数； 自由衰减振动的圆频率。设初始时刻T = 0时，初始位移 x0，初始速度为v0，则程序ex731.m%分别设 =0.1和1，wn=10，x0=1，v0=1，计算终止时间tf=2clear,wn=10;tf=2;for i=1:2 if i=1 x0=1;v0=0; else x0=0;v0=1;end %设置两组速度和位移初始条件 for j=1:10 zeta(j)=0.1*j; % 对不同的 w

46、d(j)=wn*sqrt(1-zeta(j)2); % 求wd a=sqrt(wn*x0*zeta(j)+v0)2+(x0*wd(j)2)/wd(j); % 求振幅A phi=atan2(wd(j)*x0,v0+zeta(j)*wn*x0); % 用atan2是为了求四象限相角 t=0:tf/1000:tf; % 设定自变量数组 x(j,:)=a*exp(-zeta(j)*wn*t).*sin(wd(j)*t+phi); % 求方程的解 /过渡过程 end subplot(1,2,i),plot(t,x(1,:), t,x(2,:), t,x(3,:), t,x(4,:), t,x(5,:),

47、. % 绘图 t,x(6,:), t,x(7,:), t,x(8,:), t,x(9,:), t,x(10,:) set(gcf,color,w) % 置图形背景色为白色 grid on,pauseendfigure(2),mesh(x) % 画三维图两组初始条件下 不同阻尼比对振动模态的影响三维图数据的统计描述例1：学生的身高和体重学校随机抽取100名学生，得身高和体重的数据（如下）第八章 在数理统计中的应用172 75 169 55 169 64 171 65 167 47171 62 168 67 165 52 169 62 168 65166 62 168 65 164 59 170 58 165 64160 55 175 67 175 74 172 64 168 57155 57 176 64 172 69 169 58 176 57173 58 168