-matlab在理工課程中的應(yīng)用

1、MATLAB在理工課程中的應(yīng)用第五章 在高等數(shù)學(xué)中的應(yīng)用5.1 函數(shù)、極限和導(dǎo)數(shù)一、單變量函數(shù)的計算和繪圖例5-1-1設(shè) 要求以0.01秒為間隔，求出y的151個點(diǎn)，并求出其導(dǎo)數(shù)的值和曲線。解： 建模 直接用程序文件編程的方法 編成函數(shù)文件，由主程序調(diào)用的方法t=0:.01:1.5; %設(shè)定自變量數(shù)組tw=4*sqrt(3); %固定頻率y=sqrt(3)/2*exp(-4*t).*sin(w*t+pi/3);%注意用數(shù)組運(yùn)算式subplot(2,1,1),plot(t,y),grid %繪制曲線并加上坐標(biāo)網(wǎng)格title(繪圖示例),xlabel(時間t),ylabel(y(t) % 加標(biāo)注%

2、 求導(dǎo)數(shù)并繪制導(dǎo)數(shù)曲線，注意數(shù)組求導(dǎo)數(shù)后其長度減少1Dy=diff(y);subplot(2,1,2),plot(t(1:length(t)-1),Dy),gridylabel(Dy(t) % 加標(biāo)注編成函數(shù)文件：%主程序(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 %繪制曲線并加上坐標(biāo)網(wǎng)格subplot(2,1,2),plot(t(1:length(t)-1),Dy),gridylabel(Dy(t) % 加標(biāo)注%函數(shù)文件ex511bffunc

3、tion xvalues=ex511bf(tvalues,w)xvalues=sqrt(3)/2*exp(-4*tvalues).*sin(w*tvalues+pi/3); 從本例來看，第二種方法似乎更麻煩一些，但它具有模塊化的特點(diǎn)。當(dāng)程序中要反復(fù)多次調(diào)用此函數(shù)，而且輸入不同的自變量時，利用函數(shù)文件可大大簡化編程。例5-1-2繪制極坐標(biāo)系下曲線并討論參數(shù)a,b,n的影響。為了便于比較，編制一個能分別畫兩個圖形的程序，采用for循環(huán)。 可從中看出用循環(huán)的技巧%繪制極坐標(biāo)系下曲線theta=0:0.1:2*pi; %產(chǎn)生極角向量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); %極坐標(biāo)方程 subplot(1,2,i),polar(theta,rho(i,:);%極坐標(biāo)系繪圖enda=2 b=pi/4 n=2(4葉玫瑰線)a=2 b=0 n=3(3葉玫瑰線)二、參變方程表示的函數(shù)的計算和繪圖例5-1-3擺線的繪制。當(dāng)圓輪在平面上滾動時，其圓周上任一點(diǎn)所畫的軌跡稱為擺線。 這一點(diǎn)在圓內(nèi) 內(nèi)擺線 這一點(diǎn)在圓外 外擺線rtARr=1XYOr為輪半徑，R為點(diǎn)半徑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坐標(biāo)保持等比例顯示擺線的繪制例5-1-4三次拋物線的方程為 試探討參數(shù)a和c對其圖形的影響。三、曲線族的繪制x=-2:0.1:2; %給定x數(shù)組，確定范圍及取點(diǎn)密度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坐標(biāo)等比例并確定范圍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二次曲面的方程如下 要求討論參數(shù)a,b,c對其形狀的影響，并畫出其圖形。a=input(a=);b=input(b=);c=input(c=);d=input(d=);N=input(N=); %輸入?yún)?shù)，N為網(wǎng)格線數(shù)目xgrid=linspace(-abs(a),abs(a),N); %建立x網(wǎng)格坐標(biāo)ygrid=linspace(-abs(b),abs(b),N); %建立y網(wǎng)格坐標(biāo)x,y=meshgrid(xgrid,ygrid) %確定NXN個

7、點(diǎn)的x,y網(wǎng)格坐標(biāo)z=c*sqrt(d-y.*y/b/b-x.*x/a/a);u=1; %u=1,表示z要取正負(fù)值z1=real(z); %取z的實部z1for k=2:N-1 %以下7行程序的作用是取消z中含虛數(shù)的點(diǎn) 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時加畫負(fù)半面，并加坐標(biāo)軸 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列出求空間任意曲面的交線的程序 兩空間曲面聯(lián)立起來就形成一個空間的曲線的方程 mesh語句 要顯示其交點(diǎn)，須找到各個交點(diǎn)。由于數(shù)值計算是離散點(diǎn)，難以找到完全重合的點(diǎn) 采用設(shè)置門限的方法，兩曲線差值小于設(shè)定門限，就認(rèn)為是交點(diǎn)x,y=meshgrid(-2:.1:2); %確定計算和繪圖的定義域網(wǎng)格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坐標(biāo)差小于0.1的網(wǎng)格矩陣zz=r0.*z1;yy=r0.*y;xx=r0.*x; %求這些網(wǎng)格上的坐標(biāo)值，即交線坐標(biāo)plot3(xx(r0=0),yy(r0=0),zz(r0=0),*);%畫出這些點(diǎn)colormap(gray),hold off兩曲面的交線 如果想改變曲面方程，可以改動第二行和第三行，但這樣的程序不是通用的。最好是程序運(yùn)行時能向用戶提問，允許用戶輸入曲面方程 用到字符串

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); %確定計算和繪圖的定義域網(wǎng)格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坐標(biāo)差小于0.1的網(wǎng)格矩陣zz=r0.*z1;yy=r0.*y;xx=r0.*x; %求這些網(wǎng)格上的坐標(biāo)值，即交線坐標(biāo)plot3(xx(r0=0),yy(r0=0),zz(r0=0),*);%畫出這些點(diǎn)colormap(gray),hold off5.3 數(shù)列和極限5.4 數(shù)值方法和數(shù)值積分一、任意非線性方程f(x)=0的解例5-4-1用切線法求方程的近似數(shù)值解即求任

12、意曲線y=f(x)過零點(diǎn)的問題fzero先用fplot函數(shù)得出曲線，大致了解一下曲線形狀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二、數(shù)值定積分trapz函數(shù)三、

13、多重積分例5-4-3計算二重積分 積分區(qū)域為由x=1,y=x及y=0所圍成的閉合區(qū)域xyy=x01 建模：積分區(qū)域如圖所示clear,format compactfill(0,1,1,0,0,0,1,0,y),hold %畫出積分區(qū)域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計算三重積分 積分區(qū)域為由x=1,y=x，z=xy及z=0所圍成的閉合區(qū)域 建模： 先畫出積分區(qū)域三重積分區(qū)域圖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四、微分方程和數(shù)值積分例5-4-5用數(shù)值積分法求解下列微分方程設(shè)初始時間t0=0；終止時間tf=3;初始條件y(0)=0,dy(t)/dt=0。 建模： 先將方程化成兩個一階微分方程組，設(shè)寫成矩陣形式為clf,t0=0;tf=3*pi;x0t=0;0;t,x=ode23(ex545f,t0,tf,x0t)y=x(:,1); %數(shù)值積分解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(數(shù)值積分解,輸入量,解析解)function xdot=ex545f(t,x)u=1-(t.2)/(pi2);xdot=0 1;-1 0*x+0 1*u;主程序：函數(shù)程序：5.5 線性代數(shù)矩陣的行列式det矩陣的秩rank矩陣的跡trace第六章 在普通物理中的應(yīng)用6.1 物理數(shù)據(jù)處理 例6-1-2 寫出一個程序，能把用戶輸入的長度單位在厘米，米，千米，英寸，英尺，英里，市尺，市里之間任意轉(zhuǎn)換建模： 第一步把輸入量變換為米 第二步把米變換為輸出單位 把變換常數(shù)直接表示成一個數(shù)組，選擇單位的序號也就成了數(shù)組的下標(biāo)clear all; dis

17、p( 長度單位換算程序) fprintf(長度單位: n); %選擇輸入輸出的單位fprintf( 1) 厘米 2) 米 3) 千米 4) 英寸 n);fprintf( 5) 英尺 6) 英哩 7) 市尺 8) 市里 n);InUnits = input(選擇輸入單位編號: );OutUnits = input(選擇輸出單位編號: );% 令各種單位對米的變換常數(shù)數(shù)組 ToMeter ToMeter = 0.01, 1.00, 1000.0, 0.0254, 0.3048, 1609.3, 1/3, 500 ;FrmMeter= 1./ ToMeter;%反變換常數(shù)數(shù)組FrmMeter為ToM

18、eter數(shù)組的倒數(shù)Value = input(輸入待變換的值(0為退出): ), while( Value = 0 ) ValueinM = Value*ToMeter(InUnits); %把輸入值變?yōu)槊?NewValue = ValueinM*FrmMeter(OutUnits); % 把米變?yōu)檩敵鰡挝?fprintf(變換后的值是 %g n,NewValue); % 打印變換后的值 Value = input(輸入待變換的值(0為退出): ); % 提問下個輸入值end 6.2 力學(xué)基礎(chǔ)例6-2-1 設(shè)目標(biāo)相對于設(shè)點(diǎn)的高度為yf，給定初速，試計算物體在真空中飛行的時間和距離建模：求落點(diǎn)時

19、間tf時，需要求解一個二次線性方程由解出t，它就是落地時間tf,tf有兩個解，我們得到一個有效解（大值），再求clear; y0 = 0; x0 = 0; % 初始位置vMag = input(輸入初始速度 (m/s):(書上為50) ); % 輸入初始速度vDir = input( 輸入初速方向(度):(書上為40或50) );yf = input(輸入目標(biāo)高度（米）:(書上為8) ); % 輸入目標(biāo)高度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，計算落點(diǎn)tftf=max(tf); % 去除落點(diǎn)時間tf中的庸解t=0:0.1:tf,tf; % 設(shè)定時間數(shù)組，因tf不大可能被0.1整除，必須加一個tf點(diǎn)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個電荷之間的庫侖引力。 庫侖定律是描述真空中兩個靜止點(diǎn)電荷之間相互作用的實驗定律clear all;N=input(輸入電荷數(shù)目N=:);for ic=1:N fprintf(-n 對電荷 # %g/n,ic); rc=input(輸入電荷位置x y(米):); x(ic)=rc(1); % 電荷ic的x坐標(biāo) 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結(jié)果：輸入電荷數(shù)目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 由電位的表示式計算電場并畫出等電位線和電場方向建模 如果已知空間的電位分布 則空間的電場等于電位場的負(fù)梯度其中 分別為x,y,z三個方向的單位向量MATLAB設(shè)有g(shù)radient函數(shù)，它是靠數(shù)值微分，因此空間觀測點(diǎn)應(yīng)該取得密一些，已獲得較高的精度fprintf(輸入電位分布方程 V(x,y) n);fprintf(例如: log(x.2 + y.2) n);

24、V = input(: ,s); % 讀入字符串 V(x,y)NGrid = 20; % 繪圖的網(wǎng)格線數(shù)xMax = 5; % 繪圖區(qū)從 x= -xMax 到 x= xMax yMax=5; xPlot = linspace( -xMax, xMax, NGrid); % 繪圖取的x值x,y=meshgrid(xPlot);VPlot=eval(V);ExPlot, EyPlot = gradient(-VPlot); % 電場是電位的負(fù)梯度clf; subplot(1,2,1),meshc(VPlot); % 畫含等高線的三維曲面set(gcf,color,w) % 置圖形背景色為白色xla

25、bel(x); ylabel(y); zlabel(電位);% 規(guī)定等高線圖的范圍及比例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 穩(wěn)恒磁場例6-5-1用畢奧薩伐定律計算電流環(huán)產(chǎn)生的磁場zId

26、ldBdBPrI0真空的磁導(dǎo)率0410-7H/m。任一電流元Idl在空間任一點(diǎn)P處所產(chǎn)生的磁感應(yīng)強(qiáng)度為clear all; % 清工作空間及變量初始化mu0 = 4*pi*1e-7; % 真空導(dǎo)磁率 (T*m/A)I0 = 5.0; % 環(huán)中電流(A)Rh = input(輸入環(huán)半徑Rh(m):(書上取2) );C0 = mu0/(4*pi) * I0; % 歸并常數(shù)xMax = 3; yMax = 3; % 規(guī)定圖的范圍NGx = 21; NGy = 21; % 規(guī)定觀測點(diǎn)網(wǎng)格線數(shù)x=linspace(-xMax, xMax, NGx); % 確定觀測點(diǎn)的x,y坐標(biāo)數(shù)組y=linspace(

27、-yMax, yMax, NGy);Nh = 20; % 電流環(huán)分段數(shù) % 計算每段的端點(diǎn),環(huán)在x=0平面上,其坐標(biāo)x1,x2均為零theta0 = linspace(0,2*pi, Nh+1); % 環(huán)的圓周角分段theta1 = theta0(1:Nh);y1 = Rh*cos(theta1); % 環(huán)各段向量的起點(diǎn)坐標(biāo)y1,z1z1 = Rh*sin(theta1);theta2 = theta0(2:Nh+1);y2 = Rh*cos(theta2); % 環(huán)各段向量的終點(diǎn)坐標(biāo)y2,z2z2 = Rh*sin(theta2); dlx = 0; % 計算環(huán)各段向量dl的三個分量dly

28、= y2-y1;dlz = z2-z1; xc = 0; % 計算環(huán)各段向量中點(diǎn)的三個坐標(biāo)分量yc = (y2+y1)/2;zc = (z2+z1)/2; % 循環(huán)計算各網(wǎng)格點(diǎn)上的B(x,y) 值for i=1:NGy for j=1:NGx % 對yz平面內(nèi)的電流環(huán)分段作元素群運(yùn)算,先算環(huán)上某段與觀測點(diǎn)之間的向量r rx = x(j) - xc; ry = y(i) - yc; rz = -zc; % 觀測點(diǎn)在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); % 在圖上畫出電流環(huán)xlabel(x); ylabel(y);hold off; 電流環(huán)產(chǎn)生的磁場分布圖6.6 振動與波例6-6-1振動的合成及拍

30、頻現(xiàn)象 分別輸入兩個正弦波的振幅、相位及頻率，觀察其合成結(jié)果，特別是兩個信號的頻率接近時產(chǎn)生的拍頻現(xiàn)象建模： 兩個振動相加：兩個同方向的振動分別為：用三角函數(shù)關(guān)系，可求出振動方程t=0:0.001:10;%10秒分10000個點(diǎn)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 光學(xué)例6-7-1兩點(diǎn)（雙縫）光干涉圖案 單色光通過兩個窄縫射向屏幕，相當(dāng)于位置不同的兩個同頻同相光源向屏幕照射的疊和。由于到達(dá)屏幕各點(diǎn)的距離不同，引起相位差，有的點(diǎn)加強(qiáng)，有的點(diǎn)抵消，造成干涉現(xiàn)象。yysOzL2L1S2S1光縫屏幕雙縫干涉示意圖振幅光強(qiáng)兩光源到屏幕的距離光程差相位差clcLambda = input(輸入光的波長(單位為 nm):

32、 (書上取500) );Lambda = Lambda * 1e-9; % 將nm換 變?yōu)?md = input(輸入兩個縫的間距 (單位為 mm): (書上取2) );d = d * 0.001; % 將mm 變換為 mZ = input(輸入縫到屏的距離 (單位為 m): (書上取1) );yMax = 5*Lambda*Z/d; xs = yMax; % 設(shè)定圖案的y,x 向范圍Ny=101;ys = linspace(-yMax,yMax,Ny); % y方向分成101點(diǎn)for i=1:Ny % 對屏上全部點(diǎn)進(jìn)行循環(huán)計算 % 計算第一和第二個光源到屏上各點(diǎn)的距離 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; % 計算該點(diǎn)光強(qiáng)（設(shè)兩束光強(qiáng)相同)end% 在屏上畫出圖象% clf; figure(gcf); % 清圖形窗，將它移到前面NCLevels = 255; % 確定用的灰度等級% 定標(biāo)：使最大光強(qiáng)(4.0)對應(yīng)于最大灰度級（白色）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 向的光強(qiáng)變化曲線波長(m)0.0000005光縫距離(m)0.002光柵到屏幕距離(m)1雙縫干涉條紋及光強(qiáng)分布第七章 在力學(xué)、機(jī)械中的應(yīng)用例7-1-4四連桿機(jī)構(gòu)如圖所示：輸入桿L1的轉(zhuǎn)角 求輸出桿L3轉(zhuǎn)角隨時間的變化規(guī)律，并求其角速度和角加速度。建模：四連桿機(jī)構(gòu)的運(yùn)動方程X和Y方向的長度關(guān)系確定為（1）（2）從上述兩個方程中削去 ，便可劃成一個只包括 和 的方程，給定 ，可求出滿足此方程的由（2）式得 將（1）中的 代以 得出求得 、 和 后，就可以根據(jù)桿1的角速度求出

35、桿3的角速度為了求能使f=0的 值，利用fzero函數(shù)，把f=f（ ） 定義為函數(shù)文件ex714f.m 利用global命令定義了桿件的長度參數(shù) 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(當(dāng)前角theta1=(書上取0) );w1

36、=input(桿1角速度w1=(書上取100) );theta3=input(對應(yīng)于theta1的theta3近似值=(書上取1) );th1=theta1;theta3=fzero(ex714f,theta3); % 求初始輸出theta3theta2 = asin( L3*sin(theta3)- L1*sin(theta1)/L2); % 計算三角關(guān)系w3 = L1*w1*cos(pi/2-theta1+theta2)/ (L3*cos(theta3-pi/2-theta2) ex714b.m 計算輸入輸出桿件的角位置關(guān)系和輸出角速度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(對應(yīng)于theta1最小值處的theta3近似值= (書上取1) 弧度);dt = 2*pi/180/w1; % 桿1轉(zhuǎn)2度對應(yīng)的時間增量th1=theta1(1);theta3(1)=fzero(ex714f,theta3); %求初始輸出theta3for i=2:181 th1=theta1(i); theta3(i)=fzero

38、(ex714f,theta3(i-1); % 調(diào)用ex714f函數(shù)逐次求輸出theta3endsubplot(1,2,1),plot(theta1,theta3);xlabel(theta1),ylabel(theta3),grid % 畫曲線set(gcf,color,w) % 置圖形背景色為白色w3 = diff(theta3)/dt; % 求桿3的角速度,注意求導(dǎo)數(shù)后數(shù)組長度小一subplot(1,2,2),plot(theta1(2:length(theta1),w3);grid 四連桿機(jī)構(gòu)的輸入輸出角位置關(guān)系和輸出角速度% 四連桿運(yùn)動的分析計算:求全程運(yùn)動并作動畫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(對應(yīng)于theta1最小值處的theta3近似值= (書上取1) 弧度);dt = 2*pi/180/w1; % 桿1轉(zhuǎn)2度對應(yīng)的時間增量th1=theta1(1);theta3(1)=fzero(ex714f,theta3); % 求初始輸出theta3for i=2:18

40、1 th1=theta1(i); theta3(i)=fzero(ex714f,theta3(i-1); % 調(diào)用ex714f函數(shù)逐次求輸出theta3endfigure(1)set(gcf,color,w) % 置圖形背景色為白色subplot(1,2,1),plot(theta1,theta3);xlabel(theta1),ylabel(theta3),grid %畫曲線w3 = diff(theta3)/dt; % 求桿3的角速度,注意求導(dǎo)數(shù)后數(shù)組長度小一subplot(1,2,2),plot(theta1(2:length(theta1),w3);gridpause,figure(2

41、),subplot(1,1,1),axis equal % 設(shè)定第二張圖，使縱橫比例相同axis(0,50,-20,20),axis offx0=10;y0=0; % 左支點(diǎn)坐標(biāo)x1=x0+L1*cos(theta1);y1=L1*sin(theta1); % 桿1右端點(diǎn)坐標(biāo)x2=x0+L0+L3*cos(theta3);y2=L3*sin(theta3); % 桿3左端點(diǎn)坐標(biāo)，即桿二右端點(diǎn)坐標(biāo)x3=x0+L0;y3=0; % 桿3右端點(diǎn)坐標(biāo)line(5,45,0,0) % 畫出基準(zhǔn)線h1=line(x0,x1(1),x2(1),x3(1),y0,y1(1),y2(1),y3,Linewidt

42、h,3); % 將四個端點(diǎn)相聯(lián)set(h1,erasemode,xor); % 設(shè)定更新數(shù)據(jù)的模式，把前一組數(shù)據(jù)的圖形擦除set(gcf,color,w) % 置圖形背景色為白色for i1=1:1000 i=mod(i1,180)+1; % 在前面算出的180度范圍內(nèi)的數(shù)據(jù)循環(huán)取值mod(a,b)就 是求的是a除以b的余數(shù) set(h1,Xdata,x0,x1(i),x2(i),x3,Ydata,y0,y1(i),y2(i),y3) % 更新桿端數(shù)據(jù) fft(rand(214,1); %延時， 等待一段時間，減慢更新速度 drawnow % 及時顯示end7.2材料力學(xué)例7-2-2 長為L的

43、懸臂梁，左端固定，在離固定端L1處施加力P，求它的轉(zhuǎn)角和撓度。設(shè)梁的彈性模量E=200*109N/m2,慣性矩I=2*10-4m4建模：從轉(zhuǎn)矩求轉(zhuǎn)角要經(jīng)過一次不定積分，從轉(zhuǎn)角求撓度又要經(jīng)過一次不定積分。在matlab中可通過cumsum作近似的不定積分，更為精確的函數(shù)，還可用cumtrapz函數(shù)轉(zhuǎn)角彎矩方程撓度（0 x L1）（L1x L）程序ex722.mclearL=2; P=2000; L1=1.5; % 給出常數(shù)E = 200e9; I=2e-5;x = linspace(0,L,101); dx=L/100; % 將x分100段,步長為L/100n1=L1/dx+1; % 確定x=L

44、1處對應(yīng)的下標(biāo)M1 = -P*( L1-x(1:n1); % 第一段彎矩賦值M2 = zeros(1,101-n1); % 第二段彎矩賦值（全為零）M = M1,M2; % 全梁的彎矩A = cumsum(M)*dx/(E*I); % 對彎矩積分求轉(zhuǎn)角Y = cumsum(A)*dx; % 對轉(zhuǎn)角積分求撓度subplot(3,1,1),plot(x,M),grid % 繪彎矩圖subplot(3,1,2),plot(x,A),grid % 繪轉(zhuǎn)角圖subplot(3,1,3),plot(x,Y),grid % 繪撓度圖7.3 機(jī)械振動例731 分析單自由度阻尼系統(tǒng)的阻尼系數(shù)對其固有振動模態(tài)的影

45、響建模質(zhì)量為m、阻尼系數(shù)為c、彈性系數(shù)為k的單自由度系統(tǒng)自由衰減振動時，其運(yùn)動微分方程為可改寫為式中：系統(tǒng)固有頻率； 阻尼比 小阻尼（1）時，微分方程式的解可寫為 式中：、由初始條件確定的積分常數(shù)； 自由衰減振動的圓頻率。設(shè)初始時刻T = 0時，初始位移 x0，初始速度為v0，則程序ex731.m%分別設(shè) =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 %設(shè)置兩組速度和位移初始條件 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; % 設(shè)定自變量數(shù)組 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) % 畫三維圖兩組初始條件下 不同阻尼比對振動模態(tài)的影響三維圖數(shù)據(jù)的統(tǒng)計描述例1：學(xué)生的身高和體重學(xué)校隨機(jī)抽取100名學(xué)生，得身高和體重的數(shù)據(jù)（如下）第八章 在數(shù)理統(tǒng)計中的應(yīng)用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