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

1、 MATLAB在理工課程中的應(yīng)用在理工課程中的應(yīng)用第五章第五章 在高等數(shù)學(xué)中的應(yīng)用在高等數(shù)學(xué)中的應(yīng)用5.1 5.1 函數(shù)、極限和導(dǎo)數(shù)函數(shù)、極限和導(dǎo)數(shù)一、單變量函數(shù)的計算和繪圖一、單變量函數(shù)的計算和繪圖例例5-1-1設(shè)設(shè) 要求以要求以0.01秒為間隔，求出秒為間隔，求出y的的151個點，個點，并求出其導(dǎo)數(shù)的值和曲線。并求出其導(dǎo)數(shù)的值和曲線。334234tsineyt 解：解： 建模建模 直接用程序文件編程的方法直接用程序文件編程的方法 編成函數(shù)文件，由主程序調(diào)用的方編成函數(shù)文件，由主程序調(diào)用的方法法t=0:.01:1.5; %設(shè)定自變量數(shù)組設(shè)定自變量數(shù)組tw=4*sqrt(3); %固定頻率固定

2、頻率y=sqrt(3)/2*exp(-4*t).*sin(w*t+pi/3);%注意用數(shù)組運算式注意用數(shù)組運算式subplot(2,1,1),plot(t,y),grid %繪制曲線并加上坐標(biāo)網(wǎng)格繪制曲線并加上坐標(biāo)網(wǎng)格title(繪圖示例繪圖示例),xlabel(時間時間t),ylabel(y(t) % 加標(biāo)注加標(biāo)注% 求導(dǎo)數(shù)并繪制導(dǎo)數(shù)曲線，注意數(shù)組求導(dǎo)數(shù)后其長度減求導(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)注加標(biāo)注編成函數(shù)文件：編成函數(shù)文件

3、：%主程序主程序(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)繪制曲線并加上坐標(biāo)網(wǎng)格格subplot(2,1,2),plot(t(1:length(t)-1),Dy),gridylabel(Dy(t) % 加標(biāo)注加標(biāo)注%函數(shù)文件函數(shù)文件ex511bffunction xvalues=ex511bf(tvalues,w)xvalues=sqrt(3)/2*exp(-4*tvalues).*sin(w*tvalues+pi/3)

4、;v 從本例來看，第二種方法似乎更麻從本例來看，第二種方法似乎更麻煩一些，但它具有模塊化的特點。煩一些，但它具有模塊化的特點。當(dāng)程序中要反復(fù)多次調(diào)用此函數(shù)，當(dāng)程序中要反復(fù)多次調(diào)用此函數(shù)，而且輸入不同的自變量時，利用函而且輸入不同的自變量時，利用函數(shù)文件可大大簡化編程。數(shù)文件可大大簡化編程。例例5-1-2繪制極坐標(biāo)系下曲線繪制極坐標(biāo)系下曲線并討論參數(shù)并討論參數(shù)a,b,n的影響。的影響。nbcosaq為了便于比較，編制一個能分別畫為了便于比較，編制一個能分別畫兩個圖形的程序，采用兩個圖形的程序，采用for循環(huán)。循環(huán)。q 可從中看出用循環(huán)的技巧可從中看出用循環(huán)的技巧%繪制極坐標(biāo)系下曲線繪制極坐標(biāo)系下

5、曲線theta=0:0.1:2*pi; %產(chǎn)生極角向量產(chǎn)生極角向量for i=1:2 a(i)=input(a=);b(i)=input(b=);n(i)=input(n=) rho(i,:)=a(i)*cos(b(i)+n(i)*theta); %極坐標(biāo)方程極坐標(biāo)方程 subplot(1,2,i),polar(theta,rho(i,:);%極坐標(biāo)系繪圖極坐標(biāo)系繪圖enda=2 b=pi/4 n=2(4葉玫瑰線葉玫瑰線)a=2 b=0 n=3(3葉玫瑰線葉玫瑰線)二、參變方程表示的函數(shù)的計算和繪圖二、參變方程表示的函數(shù)的計算和繪圖例例5-1-3擺線擺線的繪制。的繪制。當(dāng)圓輪在平面上滾動時，其

6、圓周上任一點當(dāng)圓輪在平面上滾動時，其圓周上任一點所畫的軌跡稱為所畫的軌跡稱為擺線擺線。 這一點在圓內(nèi)這一點在圓內(nèi) 內(nèi)擺線內(nèi)擺線 這一點在圓外這一點在圓外 外擺線外擺線rtARr=1XYOtcosRrytsinRrtxAAr為輪半徑，為輪半徑，R為點半徑為點半徑t=0:0.1:10;r=input(r=),R=input(R=)x=r*t-R*sin(t);y=r-R*cos(t);plot(x,y),hold on,axis equal %x,y坐標(biāo)保持坐標(biāo)保持等比例顯示等比例顯示擺線的繪制擺線的繪制例例5-1-4三次拋物線的方程為三次拋物線的方程為 試探討參數(shù)試探討參數(shù)a和和c對其圖形的影響

7、。對其圖形的影響。三、曲線族的繪制三、曲線族的繪制cxaxy2x=-2:0.1:2; %給定給定x數(shù)組，確定范圍及取點密度數(shù)組，確定范圍及取點密度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)等比例并坐標(biāo)等比例并確定范圍確定范圍subplot(1,2,2)for a=-3:3 plot(x,a*x.3+x),hold on,end,grid %c=1 取不同的取不同的aaxis equal ,axi

8、s(-2 2 -3 3)c和和a取不同值時的曲線族取不同值時的曲線族5.2 5.2 空間解析幾何空間解析幾何一、曲面方程一、曲面方程例例5-2-1二次曲面的方程如下二次曲面的方程如下 要求討論參數(shù)要求討論參數(shù)a,b,c對其形狀的影響，并畫對其形狀的影響，并畫出其圖形。出其圖形。dczbyax222222a=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)

9、; %建立y網(wǎng)格坐標(biāo)x,y=meshgrid(xgrid,ygrid) %確定NXN個點的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ù)的點 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時

10、加畫負(fù)半面，并加坐標(biāo)軸 axis(-abs(a),abs(a),-abs(b),abs(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列出求空間任意曲面的交線的程序列出求空間任意曲面的交線的程序v 兩空間曲面聯(lián)立起來就形成一個空間兩空間曲面聯(lián)立起來就形成一個空間的曲線的方程的曲線的方程 mesh語句語句v 要顯示其交點，須找到各個交點。由要顯示其交

11、點，須找到各個交點。由于數(shù)值計算是離散點，難以找到完全重于數(shù)值計算是離散點，難以找到完全重合的點合的點 采用設(shè)置門限的方法，兩曲采用設(shè)置門限的方法，兩曲線差值小于設(shè)定門限，就認(rèn)為是交點線差值小于設(shè)定門限，就認(rèn)為是交點2212yxzyxz322x,y=meshgrid(-2:.1:2); %確定計算和繪圖的定義域網(wǎng)格確定計算和繪圖的定義域網(wǎng)格z1=x.*x-2*y.*y; %第一個曲面方程第一個曲面方程z2=2*x-3*y; %第二個曲面方程第二個曲面方程(平面平面)mesh(x,y,z1);hold;mesh(x,y,z2); %在一個圖上同時畫兩個在一個圖上同時畫兩個圖圖r0=(abs(z1

12、-z2)=.1; %求二個曲面求二個曲面z坐標(biāo)差小于坐標(biāo)差小于0.1的網(wǎng)格的網(wǎng)格矩陣矩陣zz=r0.*z1;yy=r0.*y;xx=r0.*x; %求這些網(wǎng)格上的坐標(biāo)值，求這些網(wǎng)格上的坐標(biāo)值，即交線坐標(biāo)即交線坐標(biāo)plot3(xx(r0=0),yy(r0=0),zz(r0=0),*);%畫出這些點畫出這些點colormap(gray),hold off兩曲面的交線兩曲面的交線q 如果想改變曲面方程，可以改動第二如果想改變曲面方程，可以改動第二行和第三行，但這樣的程序不是通用的。行和第三行，但這樣的程序不是通用的。最好是程序運行時能向用戶提問，允許最好是程序運行時能向用戶提問，允許用戶輸入曲面方程

13、用戶輸入曲面方程q 用到字符串功能和eval命令 S1=input(輸入第一個方程,s);q 在原來的z1語句改為 z1=eval(s1); q 給出定義域和間隔 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)格繪圖的定義域網(wǎng)格s1=input(輸入第

14、一個方程輸入第一個方程:,s); %輸入第一個曲面方程輸入第一個曲面方程z1=eval(s1);s2=input(輸入第二個方程輸入第二個方程:,s); %輸入第二個曲面方程輸入第二個曲面方程z2=eval(s2);mesh(x,y,z1);hold;mesh(x,y,z2); %在一個圖上同時畫兩個圖在一個圖上同時畫兩個圖r0=(abs(z1-z2)=.1; %求二個曲面求二個曲面z坐標(biāo)差小于坐標(biāo)差小于0.1的網(wǎng)格矩陣的網(wǎng)格矩陣zz=r0.*z1;yy=r0.*y;xx=r0.*x; %求這些網(wǎng)格上的坐標(biāo)值，即求這些網(wǎng)格上的坐標(biāo)值，即交線坐標(biāo)交線坐標(biāo)plot3(xx(r0=0),yy(r0=

15、0),zz(r0=0),*);%畫出這些點畫出這些點colormap(gray),hold off221yxzyxz25.3 5.3 數(shù)列和極限數(shù)列和極限5.4 5.4 數(shù)值方法和數(shù)值積分?jǐn)?shù)值方法和數(shù)值積分一、任意非線性方程一、任意非線性方程f(x)=0f(x)=0的解的解例例5-4-1用切線法求方程的近似數(shù)值解用切線法求方程的近似數(shù)值解5021023xsinxxy即求任意曲線即求任意曲線y=f(x)過零點的問題過零點的問題fzero先用先用fplot函數(shù)得出曲線，大致了解一下曲線函數(shù)得出曲線，大致了解一下曲線形狀形狀fplot(x.3+10*x.2-2*sin(x)-50,-12,5)%子程

16、序子程序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ù)值定積分二、數(shù)值定積分trapz函數(shù)函數(shù)三、多重積分三、多重積分例例5-4-3計算二重積分計算二重積分dxdyyx22 積分區(qū)域為積分區(qū)域為 由由x=1,y=x及及y=0所圍成的閉所圍

17、成的閉合區(qū)域合區(qū)域xyy=x01 建模：積分區(qū)建模：積分區(qū)域如圖所示域如圖所示clear,format compactfill(0,1,1,0,0,0,1,0,y),hold %畫出積分區(qū)域畫出積分區(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=0.1s =

18、 0.3375步長步長dx=0.01s = 0.3334例例5-4-4計算三重積分計算三重積分dxdydzzxy32 積分區(qū)域為積分區(qū)域為 由由x=1,y=x，z=xy及及z=0所圍所圍成的閉合區(qū)域成的閉合區(qū)域 建模：建模： 先畫出積分區(qū)域先畫出積分區(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)=t

19、rapz(f)*dz; ends2(k)=trapz(s1)*dyends=trapz(s2)*dx步長dx=0.1 s = 4.2081步長步長dx=0.01 s = 0.0333四、微分方程和數(shù)值積分四、微分方程和數(shù)值積分例例5-4-5用數(shù)值積分法求解下列微分方程用數(shù)值積分法求解下列微分方程212tyy設(shè)初始時間設(shè)初始時間t0=0；終止時間；終止時間tf=3 ;初始條件初始條件y(0)=0,dy(t)/dt=0。 建模：建模： 先將方程化成兩個一階微分方程組，先將方程化成兩個一階微分方程組，設(shè)設(shè)21xx 21212txx212tyy寫成矩陣形式為寫成矩陣形式為211001102110011

20、0222121txtxxxxxclf,t0=0;tf=3*pi;x0t=0;0;t,x=ode23(ex545f,t0,tf,x0t)y=x(:,1); %數(shù)值積分解數(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,plot(t,y,-,t,u,+,t,y2,o)legend(數(shù)值積分解數(shù)值積分解,輸入量輸入量,解析解解析解)function xdot=ex545f(t,x)u=1-(t.2)/(pi2);xdot=0 1;-1 0*x+0 1*u;主程序

21、：主程序：函數(shù)程序：函數(shù)程序：5.5 5.5 線性代數(shù)線性代數(shù)矩陣的行列式矩陣的行列式det矩陣的秩矩陣的秩rank矩陣的跡矩陣的跡trace第六章第六章 在普通物理中的應(yīng)用在普通物理中的應(yīng)用6.1 6.1 物理數(shù)據(jù)處理物理數(shù)據(jù)處理 例例6-1-2 6-1-2 寫出一個程序，能把用戶輸入的長度寫出一個程序，能把用戶輸入的長度單位在厘米，米，千米，英寸，英尺，英里，單位在厘米，米，千米，英寸，英尺，英里，市尺，市里之間任意轉(zhuǎn)換市尺，市里之間任意轉(zhuǎn)換 建模：建模： 第一步把輸入量變換為米第一步把輸入量變換為米 第二步把米變換為輸出單位第二步把米變換為輸出單位 把變換常數(shù)直接表示成一個數(shù)組，把變換常

22、數(shù)直接表示成一個數(shù)組，選擇單位的序號也就成了數(shù)組的下標(biāo)選擇單位的序號也就成了數(shù)組的下標(biāo)clear all; disp( 長度單位換算程序長度單位換算程序) fprintf(長度單位長度單位: n); %選擇輸入輸出的單位選擇輸入輸出的單位fprintf( 1) 厘米厘米 2) 米米 3) 千米千米 4) 英寸英寸 n);fprintf( 5) 英尺英尺 6) 英哩英哩 7) 市尺市尺 8) 市里市里 n);InUnits = input(選擇輸入單位編號選擇輸入單位編號: );OutUnits = input(選擇輸出單位編號選擇輸出單位編號: );% 令各種單位對米的變換常數(shù)數(shù)組令各種單位對

23、米的變換常數(shù)數(shù)組 ToMeter ToMeter = 0.01, 1.00, 1000.0, 0.0254, 0.3048, 1609.3, 1/3, 500 ;FrmMeter= 1./ ToMeter;%反變換常數(shù)數(shù)組反變換常數(shù)數(shù)組FrmMeter為為ToMeter數(shù)組的倒數(shù)數(shù)組的倒數(shù)Value = input(輸入待變換的值輸入待變換的值(0為退出為退出): ), while( Value = 0 ) ValueinM = Value*ToMeter(InUnits); %把輸入值變?yōu)槊装演斎胫底優(yōu)槊?NewValue = ValueinM*FrmMeter(OutUnits); % 把

24、米變?yōu)檩敵鰡挝话衙鬃優(yōu)檩敵鰡挝?fprintf(變換后的值是變換后的值是 %g n,NewValue); % 打印變換后的值打印變換后的值 Value = input(輸入待變換的值輸入待變換的值(0為退出為退出): ); % 提問下個輸入值提問下個輸入值end 6.2 6.2 力學(xué)基礎(chǔ)力學(xué)基礎(chǔ)例例6-2-16-2-1 設(shè)目標(biāo)相對于設(shè)點的高度為設(shè)目標(biāo)相對于設(shè)點的高度為y yf f，給定，給定初速，試計算物體在真空中飛行的時間初速，試計算物體在真空中飛行的時間和距離和距離建模建模： 求落點時間求落點時間tf時，需要求解一個二次線性時，需要求解一個二次線性方程方程 由由 解出解出t，它就是落地時間

25、，它就是落地時間tf,tf有兩個解，有兩個解，我們得到一個有效解（大值），再求我們得到一個有效解（大值），再求f200ygt21tsinvyf00maxtcosvxclear; y0 = 0; x0 = 0; % 初始位置初始位置vMag = input(輸入初始速度輸入初始速度 (m/s):(書上為書上為50) ); % 輸入初始速度輸入初始速度vDir = input( 輸入初速方向輸入初速方向(度度):(書上為書上為40或或50) );yf = input(輸入目標(biāo)高度（米）輸入目標(biāo)高度（米）:(書上為書上為8) ); % 輸入目標(biāo)高度輸入目標(biāo)高度vx0 = vMag*cos(vDir*

26、 (pi/180); % 計算計算x,y方向的初始速度方向的初始速度vy0 = vMag*sin(vDir* (pi/180); % wy = -9.81; wx = 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; % 設(shè)定時間數(shù)組，因設(shè)定時間數(shù)組，因tf不大可能被不大可能被0.1整除，必須加一個整除，必須加一個tf點點y = y0 + vy0*t + wy*t.2/2

27、; % 計算軌跡計算軌跡x = x0 + vx0*t + wx*t.2/2;xf = max(x),plot(x,y),grid % 計算射程，畫出軌跡計算射程，畫出軌跡set(gcf,color,w) % 置圖形背景色為白色置圖形背景色為白色 例621 無阻力拋物體的飛行軌跡6.4 6.4 靜電場靜電場 例例6-4-1計算平面上計算平面上N個電荷之間的庫個電荷之間的庫侖引力。侖引力。 庫侖定律是描述真空中兩個靜止點電荷庫侖定律是描述真空中兩個靜止點電荷之間相互作用的實驗定律之間相互作用的實驗定律212212Ryyxx3021202144RqqRqqR RR RF Fo o3012214Rxx

28、qqx xF F3012214Ryyqqy yF Fclear all;N=input(輸入電荷數(shù)目輸入電荷數(shù)目N=:);for ic=1:N fprintf(-n 對電荷對電荷 # %g/n,ic); rc=input(輸入電荷位置輸入電荷位置x y(米米):); x(ic)=rc(1); % 電荷電荷ic的的x坐標(biāo)坐標(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(i

29、c)-x(jc);yij=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é)果：結(jié)果：輸入電荷數(shù)目輸入電荷數(shù)目N=:2- 對電荷對電荷 # 1/n輸入電荷位置輸入電荷位置x y(米米):1 1輸入電荷量輸入電荷量(庫侖庫侖):0.5- 對電荷

30、對電荷 # 2/n輸入電荷位置輸入電荷位置x y(米米):2 3輸入電荷量輸入電荷量(庫侖庫侖):0.7C0 = 8.9918e+009其它電荷作用在電荷其它電荷作用在電荷 # 1上的合力為上的合力為:x-分量分量:-2.81488e+008 Ny-分量分量:-5.62976e+008 N其它電荷作用在電荷其它電荷作用在電荷 # 2上的合力為上的合力為:x-分量分量:2.81488e+008 Ny-分量分量:5.62976e+008 N 例例6-4-3 由電位的表示式計算電場并畫由電位的表示式計算電場并畫出等電位線和電場方向出等電位線和電場方向建模 如果已知空間的電位分布 則空間的電場等于電位

31、場的負(fù)梯度其中 分別為x,y,z三個方向的單位向量MATLAB設(shè)有g(shù)radient函數(shù)，它是靠數(shù)值微分，因此空間觀測點應(yīng)該取得密一些，已獲得較高的精度zy,x,VV ijk kdvdvjdydvidxdvVgradientEfprintf(輸入電位分布方程輸入電位分布方程 V(x,y) n);fprintf(例如例如: log(x.2 + y.2) n);V = input(: ,s); % 讀入字符串讀入字符串 V(x,y)NGrid = 20; % 繪圖的網(wǎng)格線數(shù)繪圖的網(wǎng)格線數(shù)xMax = 5; % 繪圖區(qū)從繪圖區(qū)從 x= -xMax 到到 x= xMax yMax=5; xPlot =

32、linspace( -xMax, xMax, NGrid); % 繪圖取的繪圖取的x值值x,y=meshgrid(xPlot);VPlot=eval(V);ExPlot, EyPlot = gradient(-VPlot); % 電場是電位的負(fù)梯度電場是電位的負(fù)梯度clf; subplot(1,2,1),meshc(VPlot); % 畫含等高線的三維曲面畫含等高線的三維曲面set(gcf,color,w) % 置圖形背景色為白色置圖形背景色為白色xlabel(x); ylabel(y); zlabel(電位電位);% 規(guī)定等高線圖的范圍及比例規(guī)定等高線圖的范圍及比例subplot(1,2,2

33、), 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 穩(wěn)恒磁場穩(wěn)恒磁場例例6-5-1用畢奧薩伐定律計算電流環(huán)產(chǎn)生的用畢奧薩伐定律計算電

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

35、網(wǎng)格線數(shù) x=linspace(-xMax, xMax, NGx); % 確定觀測點的確定觀測點的x,y坐標(biāo)數(shù)組坐標(biāo)數(shù)組 y=linspace(-yMax, yMax, NGy); Nh = 20; % 電流環(huán)分段數(shù)電流環(huán)分段數(shù) % 計算每段的端點計算每段的端點,環(huán)在環(huán)在x=0平面上平面上,其坐標(biāo)其坐標(biāo)x1,x2均為零均為零 theta0 = linspace(0,2*pi, Nh+1); % 環(huán)的圓周角分段環(huán)的圓周角分段 theta1 = theta0(1:Nh); y1 = Rh*cos(theta1); % 環(huán)各段向量的起點坐標(biāo)環(huán)各段向量的起點坐標(biāo)y1,z1 z1 = Rh*sin(th

36、eta1); theta2 = theta0(2:Nh+1); y2 = Rh*cos(theta2); % 環(huán)各段向量的終點坐標(biāo)環(huán)各段向量的終點坐標(biāo)y2,z2 z2 = Rh*sin(theta2); dlx = 0; % 計算環(huán)各段向量計算環(huán)各段向量dl的三個分量的三個分量dly = y2-y1;dlz = z2-z1; xc = 0; % 計算環(huán)各段向量中點的三個坐標(biāo)分量計算環(huán)各段向量中點的三個坐標(biāo)分量yc = (y2+y1)/2;zc = (z2+z1)/2; % 循環(huán)計算各網(wǎng)格點上的循環(huán)計算各網(wǎng)格點上的B(x,y) 值值for i=1:NGy for j=1:NGx % 對對yz平面

37、內(nèi)的電流環(huán)分段作元素群運算平面內(nèi)的電流環(huán)分段作元素群運算,先算環(huán)上某段與觀測點之間先算環(huán)上某段與觀測點之間的向量的向量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 和和 y 分量分量 dlXr_y = dlz.*rx - dlx.*rz; Bx(i,j) = sum(C0*dlXr_x./r3); % 把磁場各段的把磁場各段的x

38、和和 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)在圖上畫出電流環(huán)xlabel(x); ylabel(y);hold off; 電流環(huán)產(chǎn)生的磁場分布圖電流環(huán)產(chǎn)生的磁場分布圖6.6 6.6 振動與波振動與波例例6-6-1振動的合成及拍頻現(xiàn)象振動的合成及拍頻現(xiàn)象 分別輸入兩個正弦波的振幅

39、、相位分別輸入兩個正弦波的振幅、相位及頻率，觀察其合成結(jié)果，特別是及頻率，觀察其合成結(jié)果，特別是兩個信號的頻率接近時產(chǎn)生的拍頻兩個信號的頻率接近時產(chǎn)生的拍頻現(xiàn)象現(xiàn)象 建模：建模： 1111tsinay2222tsinay兩個振動相加：兩個振動相加：兩個同方向的振動分別為：兩個同方向的振動分別為：22211121tsinatsinayyy用三角函數(shù)關(guān)系，可求出振動方程用三角函數(shù)關(guān)系，可求出振動方程t=0:0.001:10;%10秒分秒分10000個點個點a1=input(振幅振幅1);w1=input(相位相位1);a2=input(振幅振幅2);w2=input(相位相位2);y1=a1*si

40、n(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),plot(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 光學(xué)光學(xué)例例6-7-1兩點（雙縫）光兩點（雙縫）光干涉干涉圖案圖案 單色光通過兩個窄縫射向屏幕，

41、相當(dāng)單色光通過兩個窄縫射向屏幕，相當(dāng)于位置不同的兩個同頻同相光源向于位置不同的兩個同頻同相光源向屏幕照射的疊和。由于到達(dá)屏幕各屏幕照射的疊和。由于到達(dá)屏幕各點的距離不同，引起相位差，有的點的距離不同，引起相位差，有的點加強，有的點抵消，造成干涉現(xiàn)點加強，有的點抵消，造成干涉現(xiàn)象。象。yysOzL2L1S2S1光縫屏幕雙縫干涉示意圖雙縫干涉示意圖2212zdyLs2222zdyLs21LLLL2 2cos20AA2420cosBB振幅光強兩光源到屏幕的距離光程差相位差clcLambda = input(輸入光的波長輸入光的波長(單位為單位為 nm): (書上取書上取500) );Lambda =

42、 Lambda * 1e-9; % 將將nm換換 變?yōu)樽優(yōu)?md = input(輸入兩個縫的間距輸入兩個縫的間距 (單位為單位為 mm): (書上取書上取2) );d = d * 0.001; % 將將mm 變換為變換為 mZ = input(輸入縫到屏的距離輸入縫到屏的距離 (單位為單位為 m): (書上取書上取1) );yMax = 5*Lambda*Z/d; xs = yMax; % 設(shè)定圖案的設(shè)定圖案的y,x 向范圍向范圍Ny=101;ys = linspace(-yMax,yMax,Ny); % y方向分成方向分成101點點for i=1:Ny % 對屏上全部點進行循環(huán)計算對屏上全

43、部點進行循環(huán)計算 % 計算第一和第二個光源到屏上各點的距離計算第一和第二個光源到屏上各點的距離 L1 = sqrt(ys(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; % 計算該點光強（設(shè)兩束光強相同計算該點光強（設(shè)兩束光強相同)end% 在屏上畫出圖象在屏上畫出圖象% clf; figure(gcf); % 清圖形窗，將它移到前面清圖形窗，將它移到前面NCLevels = 255; % 確定用的灰

44、度等級確定用的灰度等級% 定標(biāo)：使最大光強定標(biāo)：使最大光強(4.0)對應(yīng)于最大灰度級（白色）對應(yīng)于最大灰度級（白色）Br = (B/4.0) * NCLevels;subplot(1,4,1),image(xs,ys,Br); % 畫圖象畫圖象colormap(gray(NCLevels); % 用灰度級顏色圖用灰度級顏色圖subplot(1,4,2), plot(B(:),ys) % 畫出沿畫出沿y 向的光強變化曲線向的光強變化曲線波長(m)0.0000005光縫距離(m)0.002光柵到屏幕距離(m)1雙縫干涉條紋及光強分布雙縫干涉條紋及光強分布第七章第七章 在力學(xué)、機械中的應(yīng)用在力學(xué)、機

45、械中的應(yīng)用 例例7-1-4四連桿機構(gòu)如圖所示：四連桿機構(gòu)如圖所示： 輸入桿輸入桿L1的轉(zhuǎn)角的轉(zhuǎn)角 求輸出桿求輸出桿L3轉(zhuǎn)角隨時轉(zhuǎn)角隨時間的變化規(guī)律，并求其角速度和角加速間的變化規(guī)律，并求其角速度和角加速度。度。 建模：建模： 四連桿機構(gòu)的運動方程四連桿機構(gòu)的運動方程X和和Y方向的長度方向的長度關(guān)系確定為關(guān)系確定為0LcosLcosLcosL03322110sinLsinLsinL332211（1）（2） 從上述兩個方程中削去從上述兩個方程中削去 ，便可劃成一個只包括，便可劃成一個只包括 和和 的方程，給定的方程，給定 ，可求出滿足此方程的，可求出滿足此方程的21313由（2）式得 211332

46、L/sinLsinLsin將（1）中的 代以 得出2cos22sin10LcosLLsinLsinL1LcosL,f03322113321131 求得求得 、 和和 后，就可以根據(jù)桿后，就可以根據(jù)桿1的角速的角速度求出桿度求出桿3的角速度的角速度123233211132cosL2cosL2111233222sinL2sinLL1 為了求能使為了求能使f=0的的 值，利用值，利用fzero函數(shù)，把函數(shù)，把f=f（ ） 定義為函數(shù)文件定義為函數(shù)文件ex714f.m 利用利用global命令定義了桿件的長度參數(shù)命令定義了桿件的長度參數(shù) function y=ex714f(x) global L0 L

47、1 L2 L3 th1 y=L1.*cos(th1)+L2*sqrt(1-(L3*sin(x)-L1*sin(th1).2/L2/L2)-L3*cos(x)-L0; 33 主程序主程序ex714a.m global L0 L1 L2 L3 th1 L0=20;L1=8;L2=25; L3=20; % 輸入基線及三根桿的長度輸入基線及三根桿的長度L1,L2,L3 theta1=input(當(dāng)前角當(dāng)前角theta1=(書上取書上取0) ); w1=input(桿桿1角速度角速度w1=(書上取書上取100) ); theta3=input(對應(yīng)于對應(yīng)于theta1的的theta3近似值近似值=(書上

48、取書上取1) ); th1=theta1;theta3=fzero(ex714f,theta3); % 求初始輸出求初始輸出theta3 theta2 = asin( L3*sin(theta3)- L1*sin(theta1)/L2); % 計計算三角關(guān)系算三角關(guān)系 w3 = L1*w1*cos(pi/2-theta1+theta2)/ (L3*cos(theta3-pi/2-theta2) ex714b.m 計算輸入輸出桿件的角位置關(guān)系和輸出角速度計算輸入輸出桿件的角位置關(guān)系和輸出角速度global L0 L1 L2 L3 th1 L0=20;L1=8;L2=25; L3=20;% 輸入基

49、線及三根桿的長度輸入基線及三根桿的長度L1,L2,L3w1=input(桿桿1角速度角速度w1= (書上取書上取100) 1/秒秒);theta1=linspace(0,2*pi,181); %桿桿1每圈分為每圈分為180份，間隔份，間隔2度。度。theta3=input(對應(yīng)于對應(yīng)于theta1最小值處的最小值處的theta3近似值近似值= (書上取書上取1) 弧度弧度);dt = 2*pi/180/w1; % 桿桿1轉(zhuǎn)轉(zhuǎn)2度對應(yīng)的時間增量度對應(yīng)的時間增量th1=theta1(1);theta3(1)=fzero(ex714f,theta3); %求初始輸出求初始輸出theta3for i=

50、2:181 th1=theta1(i); theta3(i)=fzero(ex714f,theta3(i-1); % 調(diào)用調(diào)用ex714f函數(shù)逐次求輸出函數(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ù)組長度小一注意求導(dǎo)數(shù)后數(shù)組長度小一subplot(1,2,2),plot(theta1(2:le

51、ngth(theta1),w3);grid 四連桿機構(gòu)的輸入輸出角位置關(guān)系和輸出角速度四連桿機構(gòu)的輸入輸出角位置關(guān)系和輸出角速度 % 四連桿運動的分析計算四連桿運動的分析計算:求全程運動并作動畫求全程運動并作動畫global L0 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;

52、 % input(對應(yīng)于對應(yīng)于theta1最小值處的最小值處的theta3近似值近似值= (書上取書上取1) 弧度弧度);dt = 2*pi/180/w1; % 桿桿1轉(zhuǎn)轉(zhuǎn)2度對應(yīng)的時間增量度對應(yīng)的時間增量th1=theta1(1);theta3(1)=fzero(ex714f,theta3); % 求初始輸出求初始輸出theta3for i=2:181 th1=theta1(i); theta3(i)=fzero(ex714f,theta3(i-1); % 調(diào)用調(diào)用ex714f函數(shù)逐次求輸出函數(shù)逐次求輸出theta3endfigure(1)set(gcf,color,w) % 置圖形背景色為

53、白色置圖形背景色為白色subplot(1,2,1),plot(theta1,theta3);xlabel(theta1),ylabel(theta3),grid %畫曲線畫曲線w3 = diff(theta3)/dt; % 求桿求桿3的角速度的角速度,注意求導(dǎo)數(shù)后數(shù)組長度小一注意求導(dǎo)數(shù)后數(shù)組長度小一subplot(1,2,2),plot(theta1(2:length(theta1),w3);gridpause,figure(2),subplot(1,1,1),axis equal % 設(shè)定第二張圖，使縱橫比例相同設(shè)定第二張圖，使縱橫比例相同axis(0,50,-20,20),axis off

54、x0=10;y0=0; % 左支點坐標(biāo)左支點坐標(biāo)x1=x0+L1*cos(theta1);y1=L1*sin(theta1); % 桿桿1右端點坐標(biāo)右端點坐標(biāo)x2=x0+L0+L3*cos(theta3);y2=L3*sin(theta3); % 桿桿3左端點坐標(biāo)，即桿二右左端點坐標(biāo)，即桿二右端點坐標(biāo)端點坐標(biāo)x3=x0+L0;y3=0; % 桿桿3右端點坐標(biāo)右端點坐標(biāo)line(5,45,0,0) % 畫出基準(zhǔn)線畫出基準(zhǔn)線h1=line(x0,x1(1),x2(1),x3(1),y0,y1(1),y2(1),y3,Linewidth,3); % 將四個端將四個端點相聯(lián)點相聯(lián)set(h1,eras

55、emode,xor); % 設(shè)定更新數(shù)據(jù)的模式，把前一組數(shù)據(jù)的圖形擦設(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)取值度范圍內(nèi)的數(shù)據(jù)循環(huán)取值mod(a,b)就就 是求的是是求的是a除以除以b的余數(shù)的余數(shù) set(h1,Xdata,x0,x1(i),x2(i),x3,Ydata,y0,y1(i),y2(i),y3) % 更新桿端數(shù)據(jù)更新桿端數(shù)據(jù) fft(rand(214,1); %延時，延時， 等待一段時間，減慢

56、更新速度等待一段時間，減慢更新速度 drawnow % 及時顯示及時顯示end7.2材料力學(xué)材料力學(xué) 例例7-2-2 長為長為L的懸臂梁，左端固定，在的懸臂梁，左端固定，在離固定端離固定端L1處施加力處施加力P，求它的轉(zhuǎn)角和，求它的轉(zhuǎn)角和撓度。設(shè)梁的彈性模量撓度。設(shè)梁的彈性模量E=200*109N/m2,慣性矩慣性矩I=2*10-4m4 建模：建模： 從轉(zhuǎn)矩求轉(zhuǎn)角要經(jīng)過一次不定積分，從從轉(zhuǎn)矩求轉(zhuǎn)角要經(jīng)過一次不定積分，從轉(zhuǎn)角求撓度又要經(jīng)過一次不定積分。在轉(zhuǎn)角求撓度又要經(jīng)過一次不定積分。在matlab中可通過中可通過cumsum作近似的不定積作近似的不定積分，更為精確的函數(shù)，還可用分，更為精確的函

57、數(shù)，還可用cumtrapz函數(shù)函數(shù)0 xLPM1x0dxM/EIAx0AdxY轉(zhuǎn)角轉(zhuǎn)角彎矩方程彎矩方程撓度撓度（0 x L1）（L1x L） 程序程序ex722.m clear L=2; P=2000; L1=1.5; % 給出常數(shù)給出常數(shù) E = 200e9; I=2e-5; x = linspace(0,L,101); dx=L/100; % 將將x分分100段段,步長步長為為L/100 n1=L1/dx+1; % 確定確定x=L1處對應(yīng)的下標(biāo)處對應(yīng)的下標(biāo) M1 = -P*( L1-x(1:n1); % 第一段彎矩賦值第一段彎矩賦值 M2 = zeros(1,101-n1); % 第二段彎

58、矩賦值（全第二段彎矩賦值（全為零）為零） M = M1,M2; % 全梁的彎矩全梁的彎矩 A = cumsum(M)*dx/(E*I); % 對彎矩積分求轉(zhuǎn)對彎矩積分求轉(zhuǎn)角角 Y = cumsum(A)*dx; % 對轉(zhuǎn)角積分求撓度對轉(zhuǎn)角積分求撓度 subplot(3,1,1),plot(x,M),grid % 繪彎矩圖繪彎矩圖 subplot(3,1,2),plot(x,A),grid % 繪轉(zhuǎn)角圖繪轉(zhuǎn)角圖 subplot(3,1,3),plot(x,Y),grid % 繪撓度圖繪撓度圖7.3 機械振動機械振動 例例731 分析單自由度阻尼系統(tǒng)的阻尼系數(shù)分析單自由度阻尼系統(tǒng)的阻尼系數(shù)對其固有

59、振動模態(tài)的影響對其固有振動模態(tài)的影響 建模建模 質(zhì)量為質(zhì)量為m、阻尼系數(shù)為、阻尼系數(shù)為c、彈性系數(shù)為、彈性系數(shù)為k的單自由度系統(tǒng)自由的單自由度系統(tǒng)自由衰減振動時，其運動微分方程為衰減振動時，其運動微分方程為0 m kxxcx 0 2xxxn 可改寫為n式中：系統(tǒng)固有頻率； 阻尼比 mkn2mkc小阻尼（ 1）時，微分方程式的解可寫為）時，微分方程式的解可寫為 )sin(dntAextAd式中：、由初始條件確定的積分常數(shù)； 自由衰減振動的圓頻率。 設(shè)初始時刻設(shè)初始時刻T = 0時，初始位移時，初始位移 x0，初始，初始速度為速度為v0，則，則2d2d0200020)x()(xvxA20n0d0)

60、(xvxarctg2nd1程序程序ex731.m%分別設(shè)分別設(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è)置兩組速度和位移初始條件設(shè)置兩組速度和位移初始條件 for j=1:10 zeta(j)=0.1*j; % 對不同的對不同的 wd(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