《计算物理(本科)》[实验例题]

1、合肥工业大学电子科学与应用物理学院实验一实验一 实验数据的统计处理法实验数据的统计处理法 一、统计直方图一、统计直方图直方图：直方图：以区间频数和该区间作为边长画成长方形，以区间频数和该区间作为边长画成长方形，这些长方形组成的图就是这些长方形组成的图就是统计直方图统计直方图。MATLAB中实现统计直方图的函数是中实现统计直方图的函数是hist(x,y)，x：向量，测量的一组数据；：向量，测量的一组数据；y：分成区间的个数，：分成区间的个数，或自定义区间的向量，缺省默认或自定义区间的向量，缺省默认10个区间。个区间。合肥工业大学电子科学与应用物理学院MATLAB中实现算术平均值的函数是中实现算术

2、平均值的函数是mean(x) x=1:100; mean(x)ans = 50.5000 x=rand(1,100); mean(x)ans = 0.5000MATLAB中实现计算测量值标准偏差的函中实现计算测量值标准偏差的函数是数是std(x)，例如：，例如： x=1:100; std(x)ans = 29.0115合肥工业大学电子科学与应用物理学院实验二实验二 实验数据的插值法实验数据的插值法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插值进行计算插值进行计算合肥工业大学电子科学与应用物理学院实验三实验三 实验数据的拟合实验数据的拟合在在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) 运行结果：运行结果：a = 0.4900 1.2501 0.8560【例】设【

4、例】设y=span1,x, x2 ,用多项式曲线拟合如表所示的数据用多项式曲线拟合如表所示的数据 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合肥工业大学电子科学与应用物理学院【例】【例】课本中课本中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运算结果：运算结果：x = 1.0000 2.0000 1.0000实验四实验四 线性代数方程组的解法线性代数方程

5、组的解法合肥工业大学电子科学与应用物理学院【例】设初值为【例】设初值为0，精度，精度10-6，用，用Jacobi迭代法求方程组解迭代法求方程组解 10 x1 - x2 =9 - x1+10 x2 - 2x3=7 -2x2+10 x3=6MALAB实现程序实现程序:a=10 -1 0; -1 10 -2;0 -2 10;b=9;7;6;x0=0; 0; 0;jacobi(a,b,x0)合肥工业大学电子科学与应用物理学院实验五实验五 数值积分及其应用数值积分及其应用解：解：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结果：结果：II= 0.9008合肥工业大学电子科学与应用物理学院2/2/202000)(141),(lldxlqyxxyxV解：如图所示建立坐标系，那么解：如图所示建立坐标系，那么(x0,y0)处的电势为处的电势为【例【例1】设长为】设长为l=0.02(m)的杆上，均匀带有电荷的杆上，均匀带有电荷q=16*10-6(库仑库仑)，可认为，可认为是线电荷，求在杆所在平是线电荷，求在杆所在平面内任意一点的电势。面内任意一点的电势。xy(x0,y0)dxo-l/2l/2r合肥工业大学电子科学与应用物理学院在在MATL

7、AB中用中用Simpson求积法对某一点求积法对某一点(x0,y0)处的处的电势进行计算。电势进行计算。1、建立、建立(x0,y0)的积分函数的积分函数%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)计算结果：计算结果：v = 1.0173e+005合肥工业大学电子科学与应用物理学院实验六实验六 常微分方程数值解法

8、及其物理中的应用常微分方程数值解法及其物理中的应用【例【例3】阻尼振动阻尼振动：求解阻尼振动方程求解阻尼振动方程022kxdtdxcdtxdm 已知质量已知质量m=10，倔，倔强系数强系数k=10，阻尼系数，阻尼系数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)合肥工业大学电子科学与应用物理学院实验七实验七 热传导方程的数值解法热传导方程的数值解法解：在解：在MATLAB中建立中建立m文件文件jswlx_9_3_1.m【例【例1】求热传导方程】求热传导方程2222yuxutu求解的范围为正方形区域：求解的范围为正方形区域：-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 。合肥工业大学电子科学与应用物理学院%求解标准热传导方程求解标准热传导方程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) 合肥工业大学电子科学与应用物理学院实验八实验八 波动方程的数值解法波动方程的数值解法解解：在：在MATLAB中建立中建立m文件文件【例【例1】求波动方程】求波动方程222222yuxutu求解的范围为正方形区域：求解的范围为正方形区域：-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 。合肥工业大学电子科学与应用物理学院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;合肥工业大学电子科学与应用物理学院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); 合肥工业大学电子科学与应用物理学院实验九实验九 蒙特卡罗方法蒙特卡罗方法例例: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)运行结果：运行结果：Fk = 16.0000 6.8284 0 1.1716 0 1.1716 0 6.8284合肥工业大学电子科学与应用物理学院【例】研究连续信号【例】研究连续信号 f(t)=