malab軟件在數(shù)值模擬中的應(yīng)用

malab軟件在數(shù)值模擬中的應(yīng)用

1b（magneter）是一個(gè)功能強(qiáng)、效率高、科學(xué)計(jì)算和應(yīng)用的嵌入式軟件包。它是美國Mathworks公司推出的一種高性能的科技計(jì)算軟件,它集數(shù)值計(jì)算、符號運(yùn)算、圖形圖像處理、編程等多種功能于一體,應(yīng)用非常廣泛。Matlab對于學(xué)習(xí)者的編程語言基礎(chǔ)要求不高,庫函數(shù)和編程語句豐富多樣且簡單易學(xué),在數(shù)據(jù)可視化上也有獨(dú)特的優(yōu)勢。學(xué)習(xí)者不需要投入太多的時(shí)間在學(xué)習(xí)編程語言知識上,可以直接利用軟件提供的豐富函數(shù),編寫較簡單的程序即可解決許多普通物理學(xué)的問題。一般而言,對從實(shí)際問題中建立的方程或公式也很難用手工計(jì)算出滿意的結(jié)果,這就需要有一種方便快捷的方法來使以上問題形象化?；诖?本文探討Matlab軟件在數(shù)值模擬方面的應(yīng)用,為利用Matlab進(jìn)行計(jì)算機(jī)模擬提供一些可供借鑒的范例。1以文字生成的方式生成文本Matlab輸入命令的方式有兩種:一種是在命令窗口中直接輸入簡單的語句,這種方式適用于命令比較簡單且處理問題沒有普遍應(yīng)用性的場合;另一種工作方式為M文件的編程工作方式。M文件是一個(gè)簡單的文本文件,語法比一般的高級語言都要簡單,交互性強(qiáng);而且可以像一般文本文件那樣在任何文本編輯器中進(jìn)行編輯、存儲(chǔ)、修改和讀取。當(dāng)使用繪圖語句時(shí),Matlab自動(dòng)打開一個(gè)圖形窗口。如果已經(jīng)有圖形窗口存在,作圖命令便會(huì)使用已存在的圖形窗口。如果使用命令figure,就會(huì)打開一個(gè)新的圖形窗口。在已有圖形上繼續(xù)作圖的指令是holdon;取消這種功能的指令是holdoff。1.1二維構(gòu)圖的對象1pool屬性二維圖形繪圖命令中最基本的指令就是plot。如果輸入兩個(gè)矢量x、y,則plot(x,y)產(chǎn)生的是y相對于x的圖形。2fpot句子當(dāng)已知函數(shù)關(guān)系的時(shí)候,采用fplot語句可以更加快速、精確地繪制出指定區(qū)間上的函數(shù)圖像。1.23d曲線設(shè)計(jì)1求取點(diǎn)坐標(biāo)其中x、y、z是同維數(shù)的矢量或矩陣。每組x、y、z構(gòu)成一個(gè)點(diǎn)的坐標(biāo),各點(diǎn)依次相連,形成一條曲線。如果是矩陣,則它們相應(yīng)的列構(gòu)成一條三維曲線的數(shù)據(jù)點(diǎn)坐標(biāo),所以用矩陣可以同時(shí)畫多條空間曲線。2生成數(shù)據(jù)網(wǎng)格的命令Matlab中,曲面是用xy平面的各個(gè)格點(diǎn)上的z坐標(biāo)來定義,相鄰點(diǎn)用直線連接。因此,為了描繪二元函數(shù)形成的空間曲面,首先要建立平面的數(shù)據(jù)網(wǎng)格。即根據(jù)x、y的定義域[x1,x2]、[y1,y2]生成平面上的坐標(biāo)格點(diǎn),每個(gè)格點(diǎn)對應(yīng)一組[x,y]。生成數(shù)據(jù)網(wǎng)格的命令是meshgrid,格式為[X,Y]=meshgrid(x,y),使用meshgrid命令,將兩個(gè)獨(dú)立的矢量x、y組成數(shù)據(jù)網(wǎng)格。mesh是三維網(wǎng)格作圖命令,mesh(x,y,z)畫出了每一個(gè)格點(diǎn)(x,y)上對應(yīng)的z值。surfl(X,Y,Z)是畫帶光照模式的三維曲面圖。該命令顯示一個(gè)帶陰影的曲面,融合了周圍的、散射的和鏡面反射的光照模式。要獲得較平滑的顏色過度,應(yīng)該使用有線性強(qiáng)度變化的色圖,如gray、copper、bone、pink等。1.3階微分方程的求解Matlab解常微分方程組的能力很強(qiáng)且方便,對于在普通物理學(xué)中遇到的大多數(shù)動(dòng)力學(xué)方程都可用命令ode45求解。Matlab只能解一階的常微分方程組,高階的常微分方程需要轉(zhuǎn)化成一階方程組才能求解。對于二階常微分方程F(¨x?˙x?x,t)=0F(x¨?x˙?x,t)=0,首先需要化成顯式形式x=f(˙x?x,t)x=f(x˙?x,t),然后令y(1)=x?y(2)=˙xy(1)=x?y(2)=x˙,則二階常微分方程化為兩個(gè)一階常微分方程組成的方程組dy(1)dt=y(2)dy(1)dt=y(2)和dy(2)dt=f(y(1),y(2)?t)dy(2)dt=f(y(1),y(2)?t),從而使問題得到解決。2bp是模擬點(diǎn)電壓的電勢和車間2.1[x,y,z]的電勢及數(shù)字圖件物理情景是Oxy平面上,在x=0、y=1處有一正電荷,x=0、y=-1處有一負(fù)電荷,求出此時(shí)的電勢分布情況。由電勢的表達(dá)式U=q4πεrU=q4πεr以及電勢疊加原理可以計(jì)算出2個(gè)點(diǎn)電荷電場的電勢分布,其中r=√(x-x0)2+(y-y0)2r=(x?x0)2+(y?y0)2?????????????????√。啟動(dòng)M文件編輯器,創(chuàng)建一個(gè)新的M文件,在編輯窗口中鍵入下列語句:[x,y]=meshgrid(-2∶0.1∶2,-2∶0.1∶2);%建立數(shù)據(jù)網(wǎng)格z=1./sqrt(x.^2+(y-1).^2+0.01)-1./sqrt(x.^2+(y+1).^2+0.01);%電勢的表達(dá)式surfl(x,y,z);%三維曲面繪圖shadinginterp;%平滑三維曲面colormap(gray);%設(shè)定顏色為灰色選定一系列x和y后,組成平面上的網(wǎng)格點(diǎn),再計(jì)算對應(yīng)每一點(diǎn)上的z值。-2∶0.1∶2,-2∶0.1∶2分別是選取橫坐標(biāo)與縱坐標(biāo)的一系列數(shù)值,meshgrid是生成數(shù)據(jù)網(wǎng)格的命令。z=1./sqrt(x.^2+(y-1).^2+0.01)-1./sqrt(x.^2+(y+1).^2+0.01)是場點(diǎn)(x,y)的電勢。當(dāng)場點(diǎn)在電荷處時(shí),會(huì)出現(xiàn)分母為零的情況,因此,在r里加了一個(gè)小量0.01,這樣既可以完成計(jì)算,又不會(huì)對結(jié)果的正確性造成太大影響。另外需要注意的是表達(dá)式中的“./”、“.^”是對數(shù)組運(yùn)算的算符,含義與數(shù)值運(yùn)算中的“/”、“^”相同,不同之處是后者只對單個(gè)數(shù)值變量進(jìn)行運(yùn)算,而前者對整個(gè)數(shù)組變量中的所有元素同時(shí)進(jìn)行運(yùn)算。surfl(x,y,z)函數(shù)繪出的三維圖形的表面是畫帶光照模式的三維曲面圖,圖形輸出如圖1(a)所示?？梢詫⑵鋼Q為mesh(x,y,z),其圖形輸出見圖1(b)。2.2電場線的繪制首先建立電場線的微分方程(二維情況)。因?yàn)殡妶鲋腥我稽c(diǎn)的電場方向都沿該點(diǎn)電場線的切線方向,所以滿足dydx=EyEx.dydx=EyEx.引入?yún)⒆兞縯得到dxEx=dyEy=t.dxEx=dyEy=t.設(shè)二點(diǎn)電荷位于(0,-1)和(0,1),二點(diǎn)電荷電量為q1和q2(均等于1),由庫倫定律和電場的疊加原理,得出微分方程˙x=dxdt=Ex=q1x[x2+(y+1)2]3/2+q2x[x2+(y-1)2]3/2?˙y=dydt=Ey=q1(y+1)[x2+(y+1)2]3/2+q2(y-1)[x2+(y-1)2]3/2.x˙=dxdt=Ex=q1x[x2+(y+1)2]3/2+q2x[x2+(y?1)2]3/2?y˙=dydt=Ey=q1(y+1)[x2+(y+1)2]3/2+q2(y?1)[x2+(y?1)2]3/2.解此方程就可以繪制出電場線。下面是微分方程的函數(shù)文件:functionydot=electricfun(t,y,flag,p1,p2);%p1,p2是參量,表示電量ydot=[p1*y(1)/(sqrt(y(1).^2+(y(2)+1).^2).^3)+p2*y(1)/(sqrt(y(1).^2+(y(2)-1).^2).^3);%dx/dt=Exp1*(y(2)+1)/(sqrt(y(1).^2+(y(2)+1).^2).^3)+p2*(y(2)-1)/(sqrt(y(1).^2+(y(2)-1).^2).^3));%dy/dt=Ey編寫好函數(shù)文件后,命名為electricfun.m,存在當(dāng)前路徑下,然后開始編寫解微分方程的主程序electric.m。p1=1;p2=1;%點(diǎn)電荷所帶電量axis([-2,2,-2,2]);%設(shè)定坐標(biāo)軸范圍-2≤x≤2,-2≤y≤2holdon;%圖形控制,不可擦除模式plot(0,1,′*r′);plot(0,-1,′*r′);%繪制兩源電荷a=(pi/24):pi/12:(2*pi-pi/24);%圓周上電場線起點(diǎn)所對應(yīng)的角度b=0.1*cos(a);c=0.1*sin(a);%電場線起點(diǎn)所對應(yīng)的相對坐標(biāo)c1=-1+c;c2=1+c;%把起點(diǎn)圓周的圓心放置在源電荷處b0=[bb];c0=[c1c2];%初始條件,所有電場線的起點(diǎn)橫、縱坐標(biāo)構(gòu)成了矢量b0和c0fori=1:48;%循環(huán)求解48次微分方程[t,y]=ode45(′electricfun′,[0:0.05:40],[b0(i),c0(i)],,p1,p2);%調(diào)用ode45求解,對應(yīng)一個(gè)初條件(起點(diǎn)),求解出一條電場線plot(y(:,1),y(:,2),′b′);%繪制出此條電場線end;%結(jié)束循環(huán),共繪制出48條電場線在程序中通過從圓周上取24個(gè)不同的角度(從π/24到2π-π/24,每隔π/12取1個(gè)角度),然后算出每個(gè)角度上起點(diǎn)的橫、縱坐標(biāo)值;[b1,c]和[b2,c]分別是以2個(gè)源點(diǎn)電荷為圓心