《計算物理(本科)》[實驗例題]

1、合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院實驗一實驗一 實驗數(shù)據(jù)的統(tǒng)計處理法實驗數(shù)據(jù)的統(tǒng)計處理法 一、統(tǒng)計直方圖一、統(tǒng)計直方圖直方圖：直方圖：以區(qū)間頻數(shù)和該區(qū)間作為邊長畫成長方形，以區(qū)間頻數(shù)和該區(qū)間作為邊長畫成長方形，這些長方形組成的圖就是這些長方形組成的圖就是統(tǒng)計直方圖統(tǒng)計直方圖。MATLAB中實現(xiàn)統(tǒng)計直方圖的函數(shù)是中實現(xiàn)統(tǒng)計直方圖的函數(shù)是hist(x,y)，x：向量，測量的一組數(shù)據(jù)；：向量，測量的一組數(shù)據(jù)；y：分成區(qū)間的個數(shù)，：分成區(qū)間的個數(shù)，或自定義區(qū)間的向量，缺省默認或自定義區(qū)間的向量，缺省默認10個區(qū)間。個區(qū)間。合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院MATLAB中實現(xiàn)算術(shù)平均值的函數(shù)是中實現(xiàn)算術(shù)

2、平均值的函數(shù)是mean(x) x=1:100; mean(x)ans = 50.5000 x=rand(1,100); mean(x)ans = 0.5000MATLAB中實現(xiàn)計算測量值標準偏差的函中實現(xiàn)計算測量值標準偏差的函數(shù)是數(shù)是std(x)，例如：，例如： x=1:100; std(x)ans = 29.0115合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院實驗二實驗二 實驗數(shù)據(jù)的插值法實驗數(shù)據(jù)的插值法x=-5:1:5; y=1./(1+x.2); x0=-5:0.1:5;y0=lagrange(x,y,x0);y1=1./(1+x0.2);plot(x0,y0,-r)hold onplot(x0,

3、y1,-b)211)(xxf【實例【實例】取】取n=10,用用Lagrange插值進行計算插值進行計算合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院實驗三實驗三 實驗數(shù)據(jù)的擬合實驗數(shù)據(jù)的擬合在在MATLAB中編寫程序中編寫程序:x=0.5 1.0 1.5 2.0 2.5 3.0;y=1.75 2.45 3.81 4.80 8.00 8.60;a=polyfit(x,y,2)x1=0.5:0.05:3.0;y1=a(3)+a(2)*x1+a(1)*x1.2;plot(x, y,*)hold on plot(x1, y1, -r) 運行結(jié)果：運行結(jié)果：a = 0.4900 1.2501 0.8560【例】設(shè)【

4、例】設(shè)y=span1,x, x2 ,用多項式曲線擬合如表所示的數(shù)據(jù)用多項式曲線擬合如表所示的數(shù)據(jù) x 0.5 1.0 1.5 2.0 2.5 3.0 y 1.75 2.45 3.81 4.80 8.00 8.60合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院【例】【例】課本中課本中P53例子：例子：2x1+3x2+x3=9x1+x2+x3=4x1-2x2-x3=-4MATLAB程序：程序：a=2.0 3.0 1.0 1.0 1.0 1.0 1.0 -2.0 -1.0;b=9 4 -4;x=ab運算結(jié)果：運算結(jié)果：x = 1.0000 2.0000 1.0000實驗四實驗四 線性代數(shù)方程組的解法線性代數(shù)方程

5、組的解法合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院【例】設(shè)初值為【例】設(shè)初值為0，精度，精度10-6，用，用Jacobi迭代法求方程組解迭代法求方程組解 10 x1 - x2 =9 - x1+10 x2 - 2x3=7 -2x2+10 x3=6MALAB實現(xiàn)程序?qū)崿F(xiàn)程序:a=10 -1 0; -1 10 -2;0 -2 10;b=9;7;6;x0=0; 0; 0;jacobi(a,b,x0)合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院實驗五實驗五 數(shù)值積分及其應(yīng)用數(shù)值積分及其應(yīng)用解：解：d=pi/1000;x=0:d:3*pi; y=exp(-0.5*x).*sin(x+pi/6);II=trapz(x,y)【

6、例】用梯形法求下列積分【例】用梯形法求下列積分305 . 0)6sin(dxxeIx結(jié)果：結(jié)果：II= 0.9008合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院2/2/202000)(141),(lldxlqyxxyxV解：如圖所示建立坐標系，那么解：如圖所示建立坐標系，那么(x0,y0)處的電勢為處的電勢為【例【例1】設(shè)長為】設(shè)長為l=0.02(m)的桿上，均勻帶有電荷的桿上，均勻帶有電荷q=16*10-6(庫侖庫侖)，可認為，可認為是線電荷，求在桿所在平是線電荷，求在桿所在平面內(nèi)任意一點的電勢。面內(nèi)任意一點的電勢。xy(x0,y0)dxo-l/2l/2r合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院在在MATL

7、AB中用中用Simpson求積法對某一點求積法對某一點(x0,y0)處的處的電勢進行計算。電勢進行計算。1、建立、建立(x0,y0)的積分函數(shù)的積分函數(shù)%jswlx_7_3_1f.mfunction v=jswlx_7_3_1f(x)x0=1;y0=1;v=16e-6/(0.02*4*pi*8.85e-12)./sqrt(x-x0).*(x-x0)+y0*y0);2、在、在MATLAB窗口中窗口中輸入輸入v=quad(jswlx_7_3_1f, -0.01,0.01,1e-10)計算結(jié)果：計算結(jié)果：v = 1.0173e+005合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院實驗六實驗六 常微分方程數(shù)值解法

8、及其物理中的應(yīng)用常微分方程數(shù)值解法及其物理中的應(yīng)用【例【例3】阻尼振動阻尼振動：求解阻尼振動方程求解阻尼振動方程022kxdtdxcdtxdm 已知質(zhì)量已知質(zhì)量m=10，倔，倔強系數(shù)強系數(shù)k=10，阻尼系數(shù)，阻尼系數(shù)c=2，初速度，初速度v0=0，初，初位置位置x0=10。解：化一階方程組解：化一階方程組 x=v , x(0)=10 v=-0.2v-x , v(0)=01、建立微分方程組、建立微分方程組m文件文件jswlx_8_4_3f.mfunction y=jswlx_8_4_3f(t,x)y=x(2);-0.2*x(2)-x(1)2、求解和繪圖文件、求解和繪圖文件jswlx_8_4_3.

9、m%Hf_1=figure; set(Hf_1,color,w)tDelt=0.05;T=70.0;x10=10;x20=0.0;x0=x10,x20;ts=0:tDelt:T;t,x=ode45(jswlx_8_4_3f,ts,x0)plot(t,x(:,1)xlabel(T),ylabel(X)合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院實驗七實驗七 熱傳導(dǎo)方程的數(shù)值解法熱傳導(dǎo)方程的數(shù)值解法解：在解：在MATLAB中建立中建立m文件文件jswlx_9_3_1.m【例【例1】求熱傳導(dǎo)方程】求熱傳導(dǎo)方程2222yuxutu求解的范圍為正方形區(qū)域：求解的范圍為正方形區(qū)域：-1x,y 1 ，t0 ；初始條件

10、：當初始條件：當x2+y2 0.4時，時，u(x,y,0)=1；其它處為；其它處為0。邊界條件：當邊界條件：當t0時，時，u(1,y,t)=0， u(x, 1,t)=0 。合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院%求解標準熱傳導(dǎo)方程求解標準熱傳導(dǎo)方程g=squareg; b=squareb1; c=1;a=0;f=0;d=1;p,e,t=initmesh(g); u0=zeros(size(p,2),1);ix=find(sqrt(p(1,:).2+p(2,:).2)0.4);u0(ix)=ones(size(ix);nframes=20;tlist=linspace(0,0.1,nframes);

11、u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);x=linspace(-1,1,31);y=x;unused,tn,a2,a3=tri2grid(p,t,u0,x,y);newplot;Mv = moviein(nframes);umax=max(max(u1);umin=min(min(u1);for j=1:nframes u=tri2grid(p,t,u1(:,j),tn,a2,a3); i= find(isnan(u); u(i)=zeros(size(i); surf(x,y,u); caxis(umin umax); colormap(cool), a

12、xis(-1 1 -1 1 0 1); Mv(:,j)=getframe; end movie(Mv,2) 合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院實驗八實驗八 波動方程的數(shù)值解法波動方程的數(shù)值解法解解：在：在MATLAB中建立中建立m文件文件【例【例1】求波動方程】求波動方程222222yuxutu求解的范圍為正方形區(qū)域：求解的范圍為正方形區(qū)域：-1x,y 1， t0 ；初始條件：當初始條件：當-1x,y1時時,第一類第一類u(x,y,0) =atancos(2x) 第二類第二類du (x,y,0) /dt= 3sin(x)expsin(y/2) ；邊界條件：當邊界條件：當t0時，時， u(-1,

13、y,t)=0， u(1,y,t)=0 。合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院g=squareg; b=squareb3; c=1;a=0;f=0;d=1;p,e,t=initmesh(squareg); x=p(1,:); y=p(2,:);u0=atan(cos(pi/2*x);ut0=3*sin(pi*x).*exp(sin(pi/2*y);n=31;tlist=linspace(0,5,n); uu=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d); delta=-1:0.1:1;uxy,tn,a2,a3=tri2grid(p,t,uu(:,1),delta

14、,delta);gp=tn;a2;a3;合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院newplot;M=moviein(n);umax=max(max(uu);umin=min(min(uu);for i=1:n,. if rem(i,10)=0,. end,. pdeplot(p,e,t,xydata,uu(:,i),zdata,uu(:,i),zstyle,continuous,. mesh,off,xygrid,on,gridparam,gp,colorbar,off);. axis(-1 1 -1 1 umin umax); caxis(umin umax);. M(:,i)=getframe;

15、. if i=n,. fprintf(donen);. end,.endnfps=5;movie(M,2,nfps); 合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院實驗九實驗九 蒙特卡羅方法蒙特卡羅方法例例:Monte Carlo方法求圓周率方法求圓周率Matlab程序程序:K=100;M=5000;p=0;for i=1:K for j=1:M X=rand(1,2); R2=(X(1,1)-0.5)2+(X(1,2)-0.5)2; if R2 fl=4 3 2 1 0 1 2 3; Fk=fft(fl)運行結(jié)果：運行結(jié)果：Fk = 16.0000 6.8284 0 1.1716 0 1.1716 0 6.8284合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院【例】研究連續(xù)信號【例】研究連續(xù)信號 f(t)=