MATLAB在普通物理中的應(yīng)用舉例

MATLAB在普通物理中的應(yīng)用舉例

6.1物理數(shù)據(jù)處理

6.2力學(xué)基礎(chǔ)

6.3分子物理學(xué)

6.4靜電場(chǎng)

6.5恒穩(wěn)磁場(chǎng)

6.6振動(dòng)與波

6.7光學(xué)

6.1物理數(shù)據(jù)處理

【例6-1-1】寫(xiě)出一個(gè)程序，要求該程序能把用戶輸入的攝氏溫度轉(zhuǎn)為華氏溫度，反之也可。

解：

◆建模兩種溫度單位之間的轉(zhuǎn)換公式為攝氏變?nèi)A氏華氏變攝氏在程序中要先考慮由用戶選擇轉(zhuǎn)換的方向，以確定選用的公式?！鬗ATLAB程序

k=input(′選擇1：攝氏變?nèi)A氏；選擇2：華氏變攝氏；鍵入數(shù)字：′);

Tin=input(′輸入待變換的溫度(允許輸入數(shù)組):′);

ifk==1Tout=Tin*9/5+32, %攝氏變?nèi)A氏

elseifk==2Tout=(Tin-32)*5/9, %華氏變攝氏

elsedisp(′未給轉(zhuǎn)換方向,轉(zhuǎn)換無(wú)效′),end

s=［′華氏′;′攝氏′］;

s1=［′轉(zhuǎn)換后的溫度為′,s(k,∶),num2str(Tout),′度′］

%注意這個(gè)語(yǔ)句的編寫(xiě)方法◆結(jié)論這是一個(gè)簡(jiǎn)單的例子，只涉及兩種單位之間的轉(zhuǎn)換，如果涉及很多單位之間的相互變換，就要多一些設(shè)計(jì)技巧，下面的例子將提供一種方法。

【例6-1-2】寫(xiě)出一個(gè)程序，要求該程序能把用戶輸入的長(zhǎng)度單位在厘米、米、千米、英寸、英尺、英里、市尺、市里之間任意轉(zhuǎn)換。

解：

◆建模這里采取的技巧分成兩步，第一步把輸入量變換為米，第二步再把米變換為輸出單位。另外，把變換常數(shù)直接表示為一個(gè)數(shù)組，選擇單位的序號(hào)也就成了數(shù)組的下標(biāo)，這樣，程序就比較簡(jiǎn)明易讀?！鬗ATLAB程序

clearall;disp(′長(zhǎng)度單位換算程序′)

fprintf(′長(zhǎng)度單位:＼n′);%選擇輸入輸出的單位

fprintf(′1)厘米2)米3)千米4)英寸＼n′); %＼n是換行符

fprintf(′5)英尺6)英里7)市尺8)市里＼n′);

InUnits=input(′選擇輸入單位:′);

OutUnits=input(′選擇輸出單位:′);

%設(shè)定各種單位對(duì)米的變換常數(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為T(mén)oMeter數(shù)組的倒數(shù) Value=input(′輸入待變換的值(0為退出):′); while(Value~=0)ValueinM=Value*ToMeter(InUnits);%把輸入值變?yōu)槊?/p>

NewValue=ValueinM*FrmMeter(OutUnits); %把米變?yōu)檩敵鰡挝?/p>

fprintf(′變換后的值是%g＼n′,NewValue); %打印變換后的值

Value=input(′輸入待變換的值(0為退出):′); %提問(wèn)下一個(gè)輸入值

end◆程序運(yùn)行結(jié)果程序在運(yùn)行后將向用戶提問(wèn)三次，若選輸入單位為米，輸出單位為英寸，待變換的值為2米，則變換后的值是78.7401英寸。

【例6-1-3】設(shè)給某元件加［1,2,3,4,5］V電壓，測(cè)得的電流為［0.2339,0.3812,0.5759,0.8153,0.9742］mA，求此元件的電阻。

解：

◆建模設(shè)直線的方程為y=a(1)x+a(2),待定的系數(shù)是a(1)、a(2)。將上述數(shù)據(jù)分別代入x、y，得:a(1)+a(2)=0.23392a(1)+a(2)=0.38123a(1)+a(2)=0.57594a(1)+a(2)=0.81535a(1)+a(2)=0.9742用矩陣形式表達(dá)這五個(gè)聯(lián)立方程,其系數(shù)矩陣設(shè)為datax和datay,有

datax*a(1)+ones(N,1)*a(2)=datay其中datax、datay都是N行數(shù)據(jù)列向量，這是N個(gè)一次代數(shù)方程構(gòu)成的方程組，含兩個(gè)未知數(shù)。方程數(shù)超過(guò)未知數(shù)個(gè)數(shù)，是一個(gè)超定方程組，寫(xiě)成A*a=B,其最小二乘解可以直接使用MATLAB的左除運(yùn)算符a=A＼B來(lái)求得。因此程序?yàn)?/p>

A=［datax,ones(N,1)］;B=datay;a=A＼B

a(2)的存在說(shuō)明此直線不通過(guò)原點(diǎn)。若要在過(guò)原點(diǎn)的曲線族中擬合，就要在原始方程中規(guī)定a(2)=0。把A中的第二列去掉，即令A(yù)=datax,a0=A＼B,如果需要用二次曲線（拋物線）來(lái)擬合數(shù)據(jù)，則結(jié)果為

A=［datax.^2,datax,ones(N,1)］;B=datay;a=A＼B用4.3節(jié)中的polyfit函數(shù)就更加簡(jiǎn)單，無(wú)需列寫(xiě)A。如果用n次曲線擬合，公式為an=polyfit(datax,datay,n),在知道系數(shù)多項(xiàng)式an后，求多項(xiàng)式值的語(yǔ)句為yi=polyval(an,xi)?！鬗ATLAB程序

clear datax=［1:5］′; datay=［0.2339,0.3812,0.5759,0.8153,0.9742］’; %原始數(shù)據(jù)

A=［datax,ones(5,1)］;B=datay;a=A＼B,r=1/a(1) %線性擬合

plot(datax,datay,′o′),holdon %繪出原始數(shù)據(jù)點(diǎn)圖

xi=0:0.1:5;yi=a(1)*xi+a(2);%設(shè)置51個(gè)取值點(diǎn)

A1=datax,a0=A1＼B,%通過(guò)原點(diǎn)的線性擬合

plot(xi,yi,xi,a0*xi,′:′)%繪圖

a2=polyfit(datax,datay,2);yi=polyval(a2,xi); %二次擬合

plot(xi,yi) %繪二次擬合曲線

holdoff圖6-1三種擬合結(jié)果◆程序運(yùn)行結(jié)果

a=0.1915

0.0217 a0=0.1974 a2=0.0049

0.1624 0.0556繪出的三種擬合曲線如圖6-1所示。6.2力學(xué)基礎(chǔ)

【例6-2-1】設(shè)目標(biāo)相對(duì)于射點(diǎn)的高度為yf，給定初速和射角,試計(jì)算物體在真空中飛行的時(shí)間和距離。

解:

◆建模無(wú)阻力拋射體的飛行是中學(xué)物理就解決了的問(wèn)題，本題的不同點(diǎn)是目標(biāo)和射點(diǎn)不在同一高度上，用MATLAB可使整個(gè)計(jì)算和繪圖過(guò)程自動(dòng)化。其好處是可快速地計(jì)算物體在不同初速和射角下的飛行時(shí)間和距離。關(guān)鍵在求落點(diǎn)時(shí)間tf時(shí)，需要解一個(gè)二次線性代數(shù)方程。由解出t，它就是落點(diǎn)時(shí)間tf。tf會(huì)有兩個(gè)解,我們只取其中一個(gè)有效解。再求xmax=v0cosθ0·tf◆MATLAB程序

clear;y0=0;x0=0; %初始位置

vMag=input(′輸入初始速度(m/s):′); %輸入初始速度的大小和方向

vDir=input(′輸入初始速度方向(°):′); yf=input(′輸入目標(biāo)高度（m）:′); %輸入目標(biāo)高度yf vx0=vMag*cos(vDir*(pi/180)); %計(jì)算x、y方向的初始速度

vy0=vMag*sin(vDir*(pi/180)); wy=-9.81;wx=0;%重力加速度(m/s^2) tf=roots(［wy/2,vy0,y0-yf］); %解二次線性代數(shù)方程，計(jì)算落點(diǎn)時(shí)間tf tf=max(tf); %去除tf兩個(gè)解中的庸解

t=0∶0.1∶tf; y=y0+vy0*t+wy*t.^2/2; %計(jì)算軌跡

x=x0+vx0*t+wx*t.^2/2; xf=max(x),plot(x,y), %計(jì)算射程，畫(huà)出軌跡在檢查曲線正確后,鍵入hold命令,把曲線保留下來(lái)，以便用同樣的初速，不同的射角，比較其軌跡和最大射程。◆程序運(yùn)行結(jié)果輸入初始速度(m/s):50輸入初始速度方向(°):40輸入目標(biāo)高度（m）:8 xf=241.0437當(dāng)初始速度方向?yàn)?0°時(shí),運(yùn)行程序得

xf=244.0677所得曲線如圖6-2所示?？梢酝ㄟ^(guò)圖形窗編輯來(lái)設(shè)置坐標(biāo)網(wǎng)格和進(jìn)行標(biāo)注。圖6-2無(wú)阻力拋射體的飛行軌跡

【例6-2-2】給定質(zhì)點(diǎn)沿x和y兩個(gè)方向的運(yùn)動(dòng)規(guī)律x(t)和y(t),求其運(yùn)動(dòng)軌跡,并計(jì)算其對(duì)原點(diǎn)的角動(dòng)量。

解:

◆建模本例要求用戶輸入運(yùn)動(dòng)規(guī)律的解析表示式，這需要用到字符串的輸入語(yǔ)句,應(yīng)當(dāng)在input語(yǔ)句中加上第二變?cè)鋝′,而運(yùn)行這個(gè)字符串要用eval命令。當(dāng)x(t)和y(t)都是周期運(yùn)動(dòng)時(shí),所得的曲線就是李薩如圖形。角動(dòng)量等于動(dòng)量與向徑的叉乘(crossproduct)。求速度需要用導(dǎo)數(shù),可用MATLAB的diff函數(shù)作近似導(dǎo)數(shù)計(jì)算。設(shè)角動(dòng)量為，質(zhì)點(diǎn)的動(dòng)量為，則在x-y平面上有L=x·mvy-y·mvx◆MATLAB程序

clear,closeall; %讀入字符串，它應(yīng)是滿足元素群運(yùn)算的語(yǔ)句

x=input(′x=′,′s′);y=input(′y=′,′s′); tf=input(′tf=′);m=1%設(shè)定質(zhì)量m，此處設(shè)為1 Ns=100;t=linspace(0,tf,Ns);dt=tf/(Ns-1); %分Ns個(gè)點(diǎn),求出時(shí)間增量dt xPlot=eval(x);yPlot=eval(y); %計(jì)算Ns個(gè)點(diǎn)的位置x(t)、y(t) %計(jì)算各點(diǎn)x(t)、y(t)的近似導(dǎo)數(shù)和角動(dòng)量，注意導(dǎo)數(shù)序列長(zhǎng)度比原函數(shù)少1 px=m*diff(xPlot)/dt;%px=Mdx/dtpy=m*diff(yPlot)/dt; %py=Mdy/dtLPlot=xPlot(1∶Ns-1).*py-yPlot(1∶Ns-1).*px; %求角動(dòng)量畫(huà)出軌跡及角動(dòng)量隨時(shí)間變化的曲線

clf;figure(gcf); %清圖形窗并把它前移

subplot(1,2,1),plot(xPlot,yPlot); %畫(huà)點(diǎn)的軌跡圖

axis(′equal′);grid %使兩軸比例相同

subplot(1,2,2),plot(t(1∶Ns-1),LPlot); %畫(huà)動(dòng)量矩隨時(shí)間的曲線◆程序運(yùn)行結(jié)果運(yùn)行此程序，輸入

x=t.*cos(t) y=t.*sin(t) tf=20后，得出圖6-3，讀者可看出其角動(dòng)量單調(diào)遞增。如果輸入

x=cos(2*t) y=sin(3*t) tf=5則得到圖6-4，其軌跡圖就是李薩如圖形。讀者可輸入其它形式的x(t)和y(t)，探討其結(jié)果。注意輸入式一定要滿足對(duì)t作元素群運(yùn)算的格式。圖6-3按方程x=tcos(t),y=tsin(t)畫(huà)出的軌跡及角動(dòng)量曲線圖6-4按方程x=cos(2t),y=cos(3t)畫(huà)出的軌跡及角動(dòng)量曲線

【例6-2-3】物體A(質(zhì)量為m1)在具有斜面的物體B(質(zhì)量為m2)上靠重力下滑,設(shè)斜面和地面均無(wú)摩擦力,求A沿斜面下滑的相對(duì)加速度a1和B的加速度a2,并求斜面和地面的支撐力N1及N2。

解:

◆建模分別畫(huà)出A和B的受力圖如圖6-5所示。對(duì)物體A列寫(xiě)動(dòng)力學(xué)方程，注意它的絕對(duì)加速度是a2與a1的合成:m1(a1cosθ-a2)=N1sinθm1a1sinθ=m1g-N1cosθ對(duì)物體B: m2a2=N1sinθ N2-N1cosθ-m2g=0圖6-5物體受力圖四個(gè)方程包含a1、a2、N1、N2等四個(gè)未知數(shù)，將含未知數(shù)的項(xiàng)移到左邊，常數(shù)項(xiàng)移到等式右邊，得到寫(xiě)成AX=B故可用MATLAB語(yǔ)句求出X=A＼B◆MATLAB程序

m1=input(′m1=′);m2=input(′m2=′); theta=input(′theta(度)=′); theta=theta*pi/180;g=9.81; A=［m1*cos(theta),-m1,-sin(theta),0;...m1*sin(theta),0,cos(theta),0;...

0,m2,-sin(theta),0;...

0,0,-cos(theta),1］; B=［0,m1*g,0,m2*g］′;X=A＼B; a1=X(1),a2=X(2),N1=X(3),N2=X(4)◆程序運(yùn)行結(jié)果輸入m1=2,m2=4,及theta=30,得到

a1=6.5400a2=1.8879 N1=15.1035N2=52.3200靜力學(xué)平衡和動(dòng)力學(xué)中求力與加速度關(guān)系的問(wèn)題,通常都可歸結(jié)為線性方程組的求解,只要方程組列寫(xiě)正確,用MATLAB的矩陣除法就可以方便而準(zhǔn)確地求解。

【例6-2-4】質(zhì)量為m的小球以速度u0正面撞擊質(zhì)量為M的靜止小球,假設(shè)碰撞是完全彈性的,即沒(méi)有能量損失,求碰撞后兩球的速度及它們與兩球質(zhì)量比K=M/m的關(guān)系。解:◆建模設(shè)碰撞后兩球速度都與u0同向，球m的速度為u，球M的速度為v，列出動(dòng)量守恒和能量守恒方程，引入質(zhì)量比K=M/m和無(wú)量綱速度ur=u/u0,vr=v/u0后，有動(dòng)量守恒

mu0=mu+Mv

(6-1)動(dòng)能守恒 (6-2)化為Kvr+ur=1

(6-3)

(6-4)由(6-3)得(6-5)代入(6-4)得(6-6)主動(dòng)球的能量損失為展開(kāi)并整理多項(xiàng)式(6-6),得這是一元二次代數(shù)方程,可用roots命令求根ur?！鬗ATLAB程序

clear K=logspace(-1,1,11); %設(shè)自變量數(shù)組K,從K=0.1~10,按等比取11個(gè)點(diǎn)

fori=1∶length(K) %對(duì)各個(gè)K循環(huán)計(jì)算

ur1=roots(［(1+1/K(i)),-2/K(i),(1/K(i)-1)］); %二次方程有兩個(gè)解

ur(i)=ur1(abs(ur1-1)>0.001);%去掉在1附近的庸解

end vr=(1-ur)./K; %用（6-5）式求vr,用元素群運(yùn)算

em=1-ur.*ur%主動(dòng)球損失的能量(相對(duì)值)［ur;vr］ %顯示輸出數(shù)據(jù)

semilogx(K,［ur,vr,em］),grid%繪圖◆程序運(yùn)行結(jié)果數(shù)字結(jié)果為(省略了幾行) Kurvrem 0.10000.81821.81820.3306 0.39810.43051.43050.8147 1.00000

1.00001.0000 2.5119

-0.43050.56950.8147 10.000

-0.81820.18180.3306繪出的曲線如圖6-6所示。圖6-6彈性碰撞后球速與K的關(guān)系可以看出,當(dāng)K>1時(shí),ur為負(fù),即當(dāng)靜止球質(zhì)量大于主動(dòng)球質(zhì)量時(shí)，主動(dòng)球?qū)⒒貜?。?dāng)K=1時(shí)，ur=0,即主動(dòng)球?qū)⑷縿?dòng)能傳給靜止球。當(dāng)K＜1時(shí)，ur為正，說(shuō)明主動(dòng)球?qū)⒗^續(xù)沿原來(lái)方向運(yùn)動(dòng)。在宏觀世界中很難找到完全彈性碰撞，而在微觀世界中，上述結(jié)果可以用來(lái)解釋康普頓效應(yīng)，即光子撞擊電子后，其散射光的波長(zhǎng)會(huì)變長(zhǎng)，而且波長(zhǎng)的增加量與散射角有關(guān)。因?yàn)楣庾拥膭?dòng)能為hυ，其中h=6.63e-34為普朗克常數(shù)，而υ是光子的角頻率，如上所述，在彈性碰撞后，光子損失了一些能量，必然表現(xiàn)為υ減小，即波長(zhǎng)增加。本例只考慮了正面碰撞，分析了不同質(zhì)量比的影響。要解釋康普頓效應(yīng)，必須考慮光子的散射是由側(cè)面碰撞產(chǎn)生的，這時(shí)可以把質(zhì)量比取定(光子的靜止質(zhì)量為零，應(yīng)考慮它的動(dòng)質(zhì)量hυ/c2,c為光速)，可以分析出其損失的能量(因而其波長(zhǎng))與散射方向相關(guān)，讀者可自行擴(kuò)展此程序進(jìn)行研究。6.3分子物理學(xué)

【例6-3-1】利用氣體分子運(yùn)動(dòng)的麥克斯韋速度分布律,求27℃下氮分子運(yùn)動(dòng)的速度分布曲線，并求速度在300～500m/s范圍內(nèi)的分子所占的比例，討論溫度T及分子量μ對(duì)速度分布曲線的影響。

解:◆建模麥克斯韋速度分布律為本例將說(shuō)明如何根據(jù)復(fù)雜數(shù)學(xué)公式繪制曲線，并研究單個(gè)參數(shù)的影響。先把麥克斯韋速度分布律列成一個(gè)子程序,以便經(jīng)常調(diào)用,并把一些常用的常數(shù)也放在其中,這樣主程序就簡(jiǎn)單了?！鬗ATLAB程序

(1)子程序:

functionf=mxwl(T,mu,v)

%mu為分子量,單位是公斤·摩爾-1(如氮為28×10－3)，即μ

%v為分子速度（可以是一個(gè)數(shù)組）

%T為氣體的絕對(duì)溫度

R=8.31;

%氣體常數(shù)

k=1.381*10^(-23);%玻爾茨曼常數(shù)

NA=6.022*10^23;%阿伏伽德羅常數(shù)

m=mu/NA;%分子質(zhì)量

%麥克斯韋分布率

f=4*pi*(m/(2*pi*k*T)).^(3/2).*exp(-m*v.^2./(2*k*T)).*v.*v;

(2)主程序: T=300;mu=28e-3;%給出T,muv=0∶1500; %給出自變量數(shù)組

y=mxwl(T,mu,v); %調(diào)用子程序

plot(v,y),holdon %畫(huà)出分布曲線

v1=300∶500; %給定速度范圍

y1=mxwl(T,mu,v1); %該范圍的分布

%畫(huà)出該曲線所圍區(qū)域

fill(［v1,500,300］,［y1,0,0］,′r′) trapz(y1) %求該范圍概率積分◆程序運(yùn)行結(jié)果執(zhí)行此程序所得的曲線及填充圖如圖6-7所示。積分結(jié)果為

ans=0.3763為了看出T和mu對(duì)曲線形狀的影響，可鍵入hold，再在程序中加上語(yǔ)句：

T=200;mu=28e-3; y=mxwl(T,mu,v);plot(v,y)%改變T，畫(huà)曲線

T=300;mu=2e-3; y=mxwl(T,mu,v);plot(v,y)%改變mu,畫(huà)曲線用圖形窗編輯加上標(biāo)注字符從圖6-7中可見(jiàn)，減小T，使分子的速度分布向低端移動(dòng)；減小分子量mu，使速度分布向高端移動(dòng)，這與物理概念是一致的。圖6-7麥克斯韋分布曲線

6.4靜電場(chǎng)

【例6-4-1】設(shè)電荷均勻分布在從z=-L到z=L，通過(guò)原點(diǎn)的線段上，其密度為q(單位為C/m)，試求出在x-y平面上的電位分布。

解：

◆建模點(diǎn)電荷產(chǎn)生的電位可表示為V=Q/4πrε0，它是一個(gè)標(biāo)量。其中r為電荷到測(cè)量點(diǎn)的距離。線電荷所產(chǎn)生的電位可用積分或疊加的方法來(lái)求。為此把線電荷分為N段，每段長(zhǎng)為dL(在MATLAB中,由于程序只認(rèn)英文字母,dL應(yīng)理解為ΔL)。每段上電荷為qdL，看做集中在中點(diǎn)的點(diǎn)電荷，它產(chǎn)生的電位為然后對(duì)全部電荷求和即可。把x-y平面分成網(wǎng)格，因?yàn)閤-y平面上的電位僅取決于離原點(diǎn)的垂直距離R,所以可以省略一維，只取R為自變量。把R從0～10m分成Nr+1個(gè)點(diǎn)，對(duì)每一點(diǎn)計(jì)算其電位?！鬗ATLAB程序

輸入線電荷半長(zhǎng)度L、分段數(shù)N、Nr及線電荷密度q的語(yǔ)句

E0=8.85e-12; %真空電介質(zhì)常數(shù)ε0 C0=1/4/pi/E0; %歸并常數(shù)

L0=linspace(-L,L,N+1); %將線電荷分N段

L1=L0(1∶N);L2=L0(2∶N+1); %確定每個(gè)線段的起點(diǎn)和終點(diǎn)

Lm=(L1+L2)/2;dL=2*L/N; %確定每個(gè)線段的中點(diǎn)坐標(biāo),Lm是N元數(shù)組

R=linspace(0,10,Nr+1); %將R分Nr+1點(diǎn) fork=1∶Nr+1 %對(duì)R的Nr+1點(diǎn)循環(huán)計(jì)算

Rk=sqrt(Lm.^2+R(k)^2); %測(cè)量點(diǎn)到各電荷段的向徑長(zhǎng)度，N元數(shù)組

Vk=C0*dL*q./Rk; %各電荷段在測(cè)量點(diǎn)處產(chǎn)生的電位，N元數(shù)組

V(k)=sum(Vk); %對(duì)各電荷段在測(cè)量點(diǎn)處產(chǎn)生的電位數(shù)組求和

end

［max(V),min(V)］ %顯示最大最小電位

plot(R,V),grid %繪圖◆程序運(yùn)行結(jié)果運(yùn)行此程序，輸入:

(1)q=1,L=5,N=50,Nr=50

(2)q=1,L=50,N=500,Nr=50所得電場(chǎng)的最大值和最小值分別為：（1）1.0e+010*［9.31994.1587］（2）1.0e+011*［1.34610.4159］沿R的電場(chǎng)分布如圖6-8所示。圖6-8線電荷產(chǎn)生的靜電位分布

【例6-4-2】由電位的表示式計(jì)算電場(chǎng),并畫(huà)出等電位線和電場(chǎng)方向。

解:

◆建模如果已知空間的電位分布V=V(x,y,z)則空間的電場(chǎng)等于電位場(chǎng)的負(fù)梯度其中，分別為x、y、z三個(gè)方向的單位向量。

MATLAB中設(shè)有g(shù)radient函數(shù)，它是靠數(shù)值微分的，因此空間觀測(cè)點(diǎn)應(yīng)取得密一些，以獲得較高的精度?！鬗ATLAB程序

V=input(′:′,′s′); %讀入字符串,例如log(x.^2+y.^2)

xMax=5;NGrid=20; %繪圖區(qū)從x=-xMax到x=xMax,網(wǎng)格線數(shù)

xPlot=linspace(-xMax,xMax,NGrid);%繪圖取x值［x,y］=meshgrid(xPlot); %x、y取同樣范圍,生成二維網(wǎng)格

VPlot=eval(V); %執(zhí)行輸入的字符串V(MATLAB語(yǔ)句)

［ExPlot,EyPlot］=gradient(-VPlot); %電場(chǎng)等于電位的負(fù)梯度

clf;subplot(1,2,1),meshc(VPlot); %畫(huà)含等高線的三維曲面

xlabel(′x′);ylabel(′y′);zlabel(′電位′); %規(guī)定等高線圖的范圍及比例

subplot(1,2,2),axis(［-xMaxxMax-xMaxxMax］); %建立第二個(gè)子圖

cs=contour(x,y,VPlot);%畫(huà)等高線,cs是等高線值

clabel(cs);holdon;%在等高線圖上加上編號(hào)

%在等高線圖上加上電場(chǎng)方向

quiver(x,y,ExPlot,EyPlot);%畫(huà)電場(chǎng)E的箭頭圖

xlabel(′x′);ylabel(′y′);holdoff;◆程序運(yùn)行結(jié)果在輸入電位方程V(x,y)=log(x.^2+y.^2)時(shí),得出圖6-9(a)所示的電位分布曲面和圖6-9(b)所示的電場(chǎng)分布向量圖。圖6-9V(x,y)=log(x.^2+y.^2)的電位三維立體圖和等位線及電場(chǎng)分布圖(a)電位三維立體圖；(b)等位線及電場(chǎng)分布圖◆MATLAB程序初始化,輸入環(huán)半徑Rh,環(huán)電流I0,真空導(dǎo)磁率mu0=4*pi*1e-7等的語(yǔ)句略

x=linspace(-3,3,20);y=x;

%確定觀測(cè)點(diǎn)的x、y坐標(biāo)數(shù)組

Nh=20; %電流環(huán)分段數(shù)

%計(jì)算每段的端點(diǎn),中點(diǎn)，電流元向量的各個(gè)分量

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);%theta1指起點(diǎn),theta2指終點(diǎn)

y2=Rh*cos(theta2);%環(huán)各段向量的終點(diǎn)坐標(biāo)y2,z2z2=Rh*sin(theta2);dlx=0;dly=y2-y1;dlz=z2-z1; %環(huán)各段向量dl的三個(gè)長(zhǎng)度分量

xc=0;yc=(y2+y1)/2;zc=(z2+z1)/2; %環(huán)各段向量中點(diǎn)的三個(gè)坐標(biāo)分量

%循環(huán)計(jì)算各網(wǎng)格點(diǎn)上的B(x,y)值

fori=1:NGy forj=1:NGx%對(duì)yz平面內(nèi)的電流環(huán)分段作元素群運(yùn)算,先算環(huán)上某段與觀測(cè)點(diǎn)之間的向量rrx=x(j)-xc;ry=y(i)-yc;rz=0-zc;%觀測(cè)點(diǎn)在z=0平面上

r3=sqrt(rx.^2+ry.^2+rz.^2).^3;%計(jì)算r^3dlXr-x=dly.*rz-dlz.*ry; %計(jì)算叉乘積dlXr的x和y分量

dlXr-y=dlz.*rx-dlx.*rz;Bx(i,j)=sum(C0*dlXr-x./r3); %把環(huán)各段產(chǎn)生的磁場(chǎng)分量累加

By(i,j)=sum(C0*dlXr-y./r3);endend%用quiver畫(huà)磁場(chǎng)向量圖

clf;quiver(x,y,Bx,By);

圖形標(biāo)注語(yǔ)句及在圖上畫(huà)出圓環(huán)位置的語(yǔ)句略

◆程序運(yùn)行結(jié)果運(yùn)行此程序所得結(jié)果如圖6-10所示,其中，A及A′處為線圈與x-y平面的兩個(gè)交點(diǎn)，A點(diǎn)電流方向垂直紙面向里，A′點(diǎn)電流方向垂直紙面向外。讀者可改變電流環(huán)的直徑來(lái)分析其影響，也可加上顯示各點(diǎn)磁場(chǎng)強(qiáng)度的語(yǔ)句來(lái)分析其強(qiáng)度的分布。

【例6-5-2】?jī)蓚€(gè)平行電流環(huán)之間截面上磁場(chǎng)分布的計(jì)算——亥姆霍茲線圈的驗(yàn)證。一對(duì)相同的共軸載流圓線圈,當(dāng)它們的間距正好等于線圈半徑時(shí),稱(chēng)之為亥姆霍茲線圈。計(jì)算表明,亥姆霍茲線圈軸線附近的磁場(chǎng)的大小十分均勻,而且都沿x方向。圖6-10電流環(huán)產(chǎn)生的磁場(chǎng)分布圖

解:

◆建模本題的計(jì)算模型與上例相同,只是把觀測(cè)區(qū)域取在兩線圈之間的小范圍內(nèi)。B線圈左邊的磁場(chǎng)就等于A線圈左邊的磁場(chǎng)，因此，A、B兩線圈在中間部分的合成磁場(chǎng)等于A線圈右邊的磁場(chǎng)與其左邊的磁場(chǎng)平移Rh后所得磁場(chǎng)的和。因此，只要觀測(cè)A線圈的左右區(qū)間x=［-Rh,Rh］內(nèi)的磁場(chǎng)就夠了，如圖6-11所示。圖6-11亥姆霍茲線圈觀測(cè)區(qū)◆MATLAB程序

clearall;

%初始化(給定環(huán)半徑、電流、圖形)

mu0=4*pi*1e-7;%真空導(dǎo)磁率(T*m/A)

I0=5.0;Rh=1; %這兩個(gè)常數(shù)不影響結(jié)果

C0=mu0/(4*pi)*I0;

%歸并常數(shù)

%下面三行輸入語(yǔ)句與上題不同,觀測(cè)范圍x?。?Rh,Rh］,即線圈的左右都取,%因?yàn)橐院笠袮線圈右邊的磁場(chǎng)和B線圈左邊的磁場(chǎng)相加

NGx=1;NGy=21; %設(shè)定觀測(cè)點(diǎn)網(wǎng)格數(shù)

x=linspace(-Rh,Rh,NGx);%設(shè)定觀測(cè)點(diǎn)范圍及數(shù)組

y=linspace(-Rh,Rh,NGy);%y也?。?Rh,Rh］

主程序段同例6-5-1(從Nh=20到最后一個(gè)end) Bax=Bx(∶,11∶21)+Bx(∶,1∶11); %x<0區(qū)域的磁場(chǎng)平移疊加到x>0區(qū)域

Bay=By(∶,11∶21)+By(∶,1∶11); subplot(1,2,1), mesh(x(11∶21),y,Bax);xlabel(′x′);ylabel(′y′); %畫(huà)出其Bx分布三維圖

subplot(1,2,2), plot(y,Bax),grid,xlabel(′y′);ylabel(′Bx′);◆運(yùn)行結(jié)果運(yùn)行結(jié)果如圖6-12所示?？梢钥闯?在亥姆霍茲線圈的兩個(gè)線圈之間的軸線附近相當(dāng)大的一個(gè)區(qū)域內(nèi),x方向的磁場(chǎng)強(qiáng)度Bx是非常均勻的。用相仿的語(yǔ)句也不難得知,該區(qū)域內(nèi)的y方向的磁場(chǎng)強(qiáng)度By近似為零，這留給讀者自行證明。為了看出這個(gè)區(qū)域的大小和形狀,可以用多種方法,這里用最簡(jiǎn)單的一種非圖形方法，鍵入圖6-12亥姆霍茲線圈軸線附近Bx按x、y的網(wǎng)格曲面和沿y向的分布圖(a)Bx按x、y的網(wǎng)格曲面;(b)Bx沿y向的分布圖

［q］=abs((Bax-Bax(11,6))/Bax(11,6))<0.02%磁場(chǎng)相對(duì)誤差＜2%的點(diǎn)

q=00000000000 00000000000 00000000000 00000000000 00010001000 00010001000 00011111000 11111111111 01111111110 01111111110 00111111100→x軸線01111111110011111111101111111111100011111000000100010000001000100000000000000000000000000000000000000000000000在圖6-12中各行和各列的間距都是0.05Rh，21行11列分別表示了兩個(gè)線圈間的xy截面,q=1處表明該點(diǎn)滿足上述磁場(chǎng)均勻性條件。應(yīng)該注意,滿足磁場(chǎng)均勻性條件的空間區(qū)域應(yīng)該是以x為軸線的回轉(zhuǎn)體。同樣鍵入p=abs((Bay)/Bax(11,6))<0.02，可找到Bay小于中心點(diǎn)Bax(11，6)的2%的區(qū)域。用format+顯示q和p能得到更緊湊的打印格式,用圖形方法能得到更美觀的表現(xiàn),例如，可用［j,k］=find(g);plot(j,k,′*′)語(yǔ)句實(shí)現(xiàn)。6.6振動(dòng)與波

【例6-6-1】振動(dòng)的合成及拍頻現(xiàn)象。分別輸入兩個(gè)正弦波的振幅、相位及頻率，觀察其合成的結(jié)果,特別是觀察當(dāng)兩個(gè)信號(hào)的頻率接近時(shí)產(chǎn)生的拍頻現(xiàn)象。

解:

◆建模兩個(gè)同方向的振動(dòng)y1=a1sin(ω1t+φ1)和y2=a2sin(ω2t+φ2)相加，得

y=y1+y2=a1sin(ω1t+φ1)+a2sin(ω2t+φ2)

用三角函數(shù)關(guān)系，可求出

當(dāng)ω1和ω2很接近時(shí)，成為一個(gè)很低的頻率，稱(chēng)為拍頻，從用MATLAB程序得到的圖形和聲音中可以很清楚地觀察拍頻現(xiàn)象?！鬗ATLAB程序

t=0:0.001:10;%給出時(shí)間軸上10s，分10000個(gè)點(diǎn)

%輸入兩組信號(hào)的振幅和頻率

a1=input(′振幅1=′);w1=input(′頻率1=′);

a2=input(′振幅2=′);w2=input(′頻率2=′);

y1=a1*sin(w1*t); %生成兩個(gè)正弦波

y2=a2*sin(w2*t);

y=y1+y2; %將兩個(gè)波疊加subplot(3,1,1),plot(t,y1),ylabel(′y1′) %畫(huà)出曲線

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),%產(chǎn)生聲音

sound(y2);pause(2),sound(y),pausesubplot(1,1,1)%繪圖復(fù)原◆程序運(yùn)行結(jié)果鍵入

a1=1.2;w1=300a2=1.8;w2=310運(yùn)行的結(jié)果如圖6-13所示。由于這兩個(gè)頻率非常接近，因此產(chǎn)生了差拍頻率。如有聲卡和音箱，也能聽(tīng)到這個(gè)拍頻。圖6-13拍頻現(xiàn)象

【例6-6-2】多普勒效應(yīng)的驗(yàn)證。設(shè)聲源從500m外以50m/s的速度向聽(tīng)者直線駛來(lái),其軌跡與聽(tīng)者的最小垂直距離為y0=20m，參看圖6-14，聲源的角頻率為1000rad/s,試求聽(tīng)者接收到的信號(hào)波形方程并生成其相應(yīng)的聲音。圖6-14聲源運(yùn)動(dòng)的幾何關(guān)系

解:

◆建模設(shè)聲源發(fā)出的信號(hào)為f(t)，傳到聽(tīng)者處，被聽(tīng)者接收的信號(hào)經(jīng)歷了聲音傳播的延遲，延遲時(shí)間為其中，c為音速，r為聲源與聽(tīng)者之間的距離，故接收的信號(hào)形式為(不考慮聲波的傳輸衰減)只要給出f(t)及r隨t變化的關(guān)系，即可求得f1(t)并將它恢復(fù)為聲音信號(hào)。◆MATLAB程序

x0=500;v=60;y0=30; %設(shè)定聲源運(yùn)動(dòng)參數(shù)

c=330;w=1000; %音速和頻率

t=0:0.001:30; %設(shè)定時(shí)間數(shù)組

r=sqrt((x0-v*t).^2+y0.^2);%計(jì)算聲源與聽(tīng)者的距離

t1=t-r/c;

%經(jīng)距離遲延后的等效時(shí)間

u=sin(w*t)+sin(1.1*w*t); %聲源發(fā)出的信號(hào)設(shè)為兩種頻率的合成

u1=sin(w*t1)+sin(1.1*w*t1); %聽(tīng)者接收到的信號(hào)

sound(u);pause(5);sound(u1); %將原信號(hào)和接收到的信號(hào)恢復(fù)為聲音◆程序運(yùn)行結(jié)果打開(kāi)計(jì)算機(jī)的聲音系統(tǒng),運(yùn)行此程序?qū)?huì)聽(tīng)到類(lèi)似于火車(chē)汽笛的聲音。第一聲是火車(chē)靜止時(shí)的汽笛聲,第二聲是本題中靜止的聽(tīng)者聽(tīng)到的運(yùn)動(dòng)火車(chē)的汽笛聲,它的頻率先高于原來(lái)的汽笛聲,后低于原來(lái)的汽笛聲。程序中兩個(gè)sound語(yǔ)句之間的pause(暫停)語(yǔ)句是不可少的,而且暫停的時(shí)間要足夠長(zhǎng)，以便再打開(kāi)聲音系統(tǒng)，這個(gè)量與計(jì)算機(jī)硬件有關(guān)，在作者的計(jì)算機(jī)上，這個(gè)數(shù)至少要取5。6.7光學(xué)

【例6-7-1】?jī)牲c(diǎn)（雙縫）光干涉圖案。單色光通過(guò)兩個(gè)窄縫射向屏幕，相當(dāng)于位置不同的兩個(gè)同頻同相光源向屏幕照射的疊合，由于到達(dá)屏幕各點(diǎn)的距離（光程）不同引起相位差，如圖6-15所示，疊合的結(jié)果是在有的點(diǎn)加強(qiáng)，在有的點(diǎn)抵消，造成干涉現(xiàn)象。純粹的單色光不易獲得，通常都有一定的光譜寬度，這種光的非單色性對(duì)光的干涉會(huì)產(chǎn)生何種效應(yīng)，要求用MATLAB計(jì)算并仿真這一問(wèn)題。圖6-15雙縫干涉的示意圖

解:

◆建?？紤]兩個(gè)離中心點(diǎn)距離各為d/2的相干光源S1和S2到屏幕上任意點(diǎn)的距離差引起的相位差,先分析光程差：則光程差為ΔL=L1-L2將ΔL除以波長(zhǎng)λ，并乘以2π，得到相位差。設(shè)兩束相干光在屏幕上產(chǎn)生的幅度相同，均為A0，則夾角為的兩個(gè)向量合成向量的幅度為光強(qiáng)B正比于振幅A的平方，令B=kA2,

,故有根據(jù)這些關(guān)系式，可以編寫(xiě)出計(jì)算屏幕上各點(diǎn)光強(qiáng)的程序，本書(shū)中為exn671?？紤]到光的非單色性對(duì)干涉條紋的影響，將使問(wèn)題更為復(fù)雜，此時(shí)波長(zhǎng)將不是常數(shù)，因此必須對(duì)不同波長(zhǎng)的光作分類(lèi)處理再疊加起來(lái)。假定光源的光譜寬度為中心波長(zhǎng)的±10%，并且在該區(qū)域內(nèi)均勻分布。在(0.9~1.1)λ之間，按均等間距近似取11根譜線，其波長(zhǎng)分別為相位差的計(jì)算式求出的將是11根不同譜線的11個(gè)不同相位。計(jì)算光強(qiáng)時(shí)應(yīng)把這11根譜線產(chǎn)生的光強(qiáng)疊加并取平均值，即

◆MATLAB程序(在本書(shū)程序集中將分成兩個(gè)程序exn671和exn671a)

輸入波長(zhǎng)Lambda=500nm,光縫距離d=2mm,光柵到屏幕距離z=1m yMax=5*Lambda*Z/d;xs=yMax;%設(shè)定圖案的y、x向范圍

Ny=101;ys=linspace(-yMax,yMax,Ny);%y方向分成101點(diǎn)

fori=1∶Ny%對(duì)屏上全部點(diǎn)進(jìn)行循環(huán)計(jì)算

%計(jì)算第一個(gè)和第二個(gè)光源到屏上各點(diǎn)的距離

L1=sqrt((ys(i)-d/2).^2+z^2);L2=sqrt((ys(i)+d/2).^2+z^2);Phi=2*pi*(L2-L1)/Lambda;%從距離差計(jì)算相位差

B(i,∶)=4*cos(Phi/2).^2;%計(jì)算該點(diǎn)光強(qiáng)（設(shè)兩束光強(qiáng)相同)%若考慮光譜的非單色性,把前兩句改為后四句(可在前兩句前加“%”號(hào)，而取消后四句前的“%”號(hào))，由此構(gòu)成程序exn671a

%Nl=11;dL=linspace(-0.1,0.1,Nl); %設(shè)光譜相對(duì)寬度±10%%Lambda1=Lambda*(1+dL′); %分11根譜線,波長(zhǎng)為一個(gè)數(shù)組

%Phi1=2*pi*(L2-L1)./Lambda1; %從距離差計(jì)算各波長(zhǎng)的相位差

%B(i,∶)=sum(4*cos(Phi1/