Matlab在物理上的应用举例

1、单列波 %单列波 t=0:0.001:10;A=input( 振幅A=);w=input( 频率w=);a=input(a=);y=A.*sin(w.*t+a);plot(t,y);pause(1),sound(y);ylabel(y),xlabel(t)%光的单缝衍射现象Lambda=500e-9; % a=input(a=);% 可取 0.2e-3,1e-3,2e-3 三种情况z=1%ymax=3*Lambda*z/a;%Ny=51; % ys=linspace(-ymax,ymax,Ny); % NPoints=51; % yPoint=linspace(-a/2,a/2,NPoints

2、); % for j=1:Ny %L=sqrt(ys(j)-yPoint).A2+zA2); %Phi=2*pi.*(L-z)./Lambda; %SumCos=sum(cos(Phi); %SumSin=sum(sin(Phi); %B(j)=(SumCosA2+SumSinA2)/NPointsA2; % end clf,plot(ys,B,*,ys,B);grid; % 用毕奥沙伐尔定律计算电流环产生的磁场 mu0=4*pi*1e-7; I0=5.0;Rh=1;C0=mu0/(4*pi)*I0;NGx=21;NGy=21;x=linspace(-Rh,Rh,NGx);y=linspace

3、(-3,3,20);y=x;Nh=20;theta0=linspace(0,2*pi,Nh+1);theta1=theta0(1:Nh);y1=Rh*cos(theta1);z1=Rh*sin(theta1);theta2=theta0(2:Nh+1);y2=Rh*cos(theta2);z2=Rh*sin(theta2);dlx=0;dly=y2-y1;dlz=z2-z1;xc=0;yc=(y2+y1)/2;zc=(z2+z1)/2;for i=1:NGyfor j=1:NGxrx=x(j)-xc;ry=y(i)-yc;rz=0-zc;r3=sqrt(rx.A2+ry.A2+rz.A2).A

4、3;dlXr_x=dly.*rz-dlz.*ry;dlXr_y=dly.*rx-dlx.*rz;Bx(i,j)=sum(C0*dlXr_x./r3);By(i,j)=sum(C0*dlXr_y./r3);end end clf;quiver(x,y,Bx,By);%多普勒效应x0=500;v=50;y0=20;c=330;w=1000;t=0:0.001:30;r=sqrt(x0-v*t).A2+y0.A2);t1=t-r/c;u=sin(w*t)+sin(1.1*w*t);u1=sin(w*t1)+sin(1.1*w*t1);sound(u);pause(5);sound(u1);亥姆霍兹线

5、圈clear allmu0=4*pi*1e-7;I0=5.0;Rh=1;C0=mu0/(4*pi)*I0;NGx=21;NGy=21;x=linspace(-Rh,Rh,NGx);y=linspace(-Rh,Rh,NGy);Nh=20;theta0=linspace(0,2*pi,Nh+1);theta1=theta0(1:Nh);y1=Rh*cos(theta1);z1=Rh*sin(theta1);theta2=theta0(2:Nh+1);y2=Rh*cos(theta2);z2=Rh*sin(theta2);dlx=0;dly=y2-y1;dlz=z2-z1;xc=0;yc=(y2+

6、y1)/2;zc=(z2+z1)/2;for i=1:NGyfor j=1:NGxrx=x(j)-xc;ry=y(i)-yc;rz=0-zc;r3=sqrt(rx.A2+ry.A2+rz.A2).A3;dlXr_x=dly.*rz-dlz.*ry;dlXr_y=dly.*rx-dlx.*rz;Bx(i,j)=sum(C0*dlXr_x./r3);By(i,j)=sum(C0*dlXr_y./r3);endendBax=Bx(:,11:21)+Bx(:,1:11);Bay=By(:,11:21)+By(:,1:11);subplot(1,2,1)mesh(x(11:21),y,Bax);xlab

7、el(x);ylabel(y); subplot(1,2,2), plot(y,Bax),grid,xlabel(y);ylabel(Bx);.库仑引力clear all;N=input( 电荷数目 N ： );for ic=1:N %fprintf(n 对电荷 #%gn,ic);rc=input( 电荷位置 x y (米) ： );x(ic)=rc(1);%y(ic)=rc(2);%q(ic)=input( 输入电荷量(库仑) );endE0=8.85e-12; %C0=1/(4*pi*E0);%for ic=1:NFx=0.0; Fy=0:0;%for jc=1:N %if(ic=jc)

8、%xij=x(ic)-x(jc); yij=y(ic)-y(jc);Rij=sqrt(xijA2+yijA2); %Fx=Fx+C0*q(ic)*xij/RijA3;Fy=Fy+C0*q(ic)*yij/RijA3;end end fprintf( 其它电荷作用在电荷 #%g 上的合力为： n,ic);%fprintf(x 分量：%gNn,Fx);fprintf(y 分量：%gNn,Fy)end 7 .李萨如图形% lisaru.msyms t a1 a2 w1 w2x=cos(w1.*t+a1);y=sin(w2.*t+a2);a1=input(a1=);a2=input(a2=);w1=i

9、nput(w1=);w2=input(w2=);tf=10;Ns=1000;t=linspace(0,tf,Ns);dt=tf/(Ns-1);%分 Ns 个点,求出时间增量dtxplot=eval(x);yplot=eval(y);%计算 Ns 个点的位置x(t),y(t)figure(gcf); subplot(1,2,1),for i=1:750plot(yplot(1:i),xplot(1:i); % 画点的轨迹图axis(equal); grid ;%使两轴比例相同pause(0.01)end%?%?%?u(:,s)? %?%?%?%?%?%?m1=2;m2=2;K1=16;K2=4;

10、x0=1;0;xd0=6;6;tf=10;M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2 ;u,L=eig(K,M) ;t=linspace(0,tf,101);x=zeros(2,101);for s=1:2alfa=sqrt(u(:,s)*M*u(:,s) ;u(:,s)= u(:,s)/alfa;w(j)=sqrt(L(j,j);xt=u(:,j)*(u(:,j)*M*x0*cos(w(j)*t)+u(:,j)*M*xd0/w(j)*sin(w(j)*t);x=x+xt;%?endfor r=1:2%?x1,x1?subplot(2,1,r)plot(t.x(r,:)

11、,grid;xlabel(xxx);ylabel(yyy,num2str(r);end9 .拍频% t=0:0.001:10;a1=input(?1=);w1=input(?1=);a2=input(?2=);w2=input(?2=);y1=a1*sin(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(5),

12、sound(y2);pause(5),sound(y),pause subplot(1,1,1)1 0。物块下滑 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=AB;a1=X(1), a2=X(2), N1=X(3),N2=X(4)1 1。循环

13、%xunhuanR=8.31;gama=1.4; %(注：由于在 MATLAB 中无拉丁文，所以用gama代替？)nMoles=0.5;P(1)=1e5;V(1)=0.012;WTotal=0;QTotal=0;iPoint=1;NCurve=100;PPlot=P(1);VPlot=V(1);%变量和图形初始化，输入气体的摩尔数nMoles,初始压力P (1),初始容积V (1),%气体常数R= 8.314;给定起始总功WTotal=0;点序号iPoint=1;画等温线用%的点数 NCurve=100;P V 图第一点坐标 PPlot=P(1);VPlot=V(1)T(1)=P(1)*V(1

14、)/(nMoles*R);%算出初始温度%为了进入循环，先要设两个不相等的PathType和QuitType值QuitType=5;PathType=0;while(PathType=QuitType) %在菜单上选择退出之前不断循环，% 选择路径类型或退出iPoint=iPoint+1;%下一点fprintf( 对过程gn,iPoint-1);PathType=menu(sprintf( 过程g ：选择下一路径,iPoint-1),.等压,等容 ,等温 ,绝热 ,退出);% 图形界面菜单生成语句switch PathType case 1%等压路径V(iPoint)=input( 输入新容积

15、： );P(iPoint)=P(iPoint-1);%压力不变T(iPoint)=P(iPoint)*V(iPoint)/(nMoles*R); %按新容积算出温度W=P(iPoint)*(V(iPoint)-V(iPoint-1); %计算等压过程所做的功Q=(gama*nMoles*R/(gama-1)*(T(iPoint)-T(iPoint-1);VPlot=VPlot,V(iPoint); %加上新的容积和压力点，用以绘图PPlot=PPlot,P(iPoint);case 2%等容路径P(iPoint)=input( 输入新压力： );V(iPoint)=V(iPoint-1);%容

16、积不变T(iPoint)=P(iPoint)*V(iPoint)/(nMoles*R); %按新压力算出温度W=0;%等容路径上所做的功为零Q=(nMoles*R/(gama-1)*(T(iPoint)-T(iPoint-1);VPlot=VPlot,V(iPoint);%加上绘图用的新容积和压力点PPlot=PPlot,P(iPoint); case 3%等温路径V(iPoint)=input( 输入新容积： );T(iPoint)=T(iPoint-1);%温度不变P(iPoint)=nMoles*R*T(iPoint)/V(iPoint); %按新容积求新压力W=nMoles*R*T(i

17、Point)*log(V(iPoint)/V(iPoint-1); %求所做的功 Q=W; %用元素群运算求等温路径上的P和V ,加进绘图数据中VNew=linspace(V(iPoint-1),V(iPoint),NCurve);PNew=nMoles*R*T(iPoint)./VNew;VPlot=VPlot,VNew;%将新的 V ， P 点加入绘图数据中PPlot=PPlot,PNew;case 4% 绝热路径V(iPoint)=input( 输入新容积:);P(iPoint尸P(iPoint-1)*V(iPoint-1Fgama/V(iPoint)Agama; %按新容积求新压力T(

18、iPoint尸T(iPoint-1)*V(iPoint-1)A(gama-1)/V(iPoint)A(gama-1); %按新容积求新温度W=-(nMoles*R/(gama-1)*(T(iPoint)-T(iPoint-1); % 求所做的功Q=0;VNew=linspace(V(iPoint-1),V(iPoint),NCurve);PNew=P(iPoint-1)*V(iPoint-14gama./VNew.Agama;VPlot=VPlot,VNew;%将新的 V ， P 点加入绘图数据中PPlot=PPlot,PNew; otherwise end %画出到目前为止的 PV 图 if

19、(PathType=QuitType)WTotal=WTotal+W;%将新做的功加进总功if(Q0)QTotal=QTotal+Q;endfigure(gcf);plot(V,P,o,VPlot,PPlot,-)% 图形窗移前，绘图size=axis;axis(size(1)*0.9,size(2)*1.1,size(3)*0.9,size(4)*1.1);%上两句用于美观图象，使曲线不要紧靠边框%标注语句略endendWTotalQTotaly=WTotal/QTotal1 2。振子和拍m1=2;m2=2;K1=8;K2=4;%输入各原始参数x0=1;0;xd0=6;0;tf=10;%初使

20、条件M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2 ;%构成参数矩阵u(:,s) *u,L=eig(K,M) ;%求广义特征向量和特征值t=linspace(0,tf,101);x=zeros(2,101);%时间分割和输出变量初始化for s=1:2%分别处理两特征值for j=1:2 alfa=sqrt(u(:,s)*M*u(:,s) ; u(:,s)= u(:,s)/alfa;m=sqrt(L);w(j)=m(j,j)% 分别求对应于两特征值的分量% 解耦后的向量 %把特征向量归一化xt=u(:,j)*(u(:,j)*M*x0*cos(w(j)*t)+u(:,j)*M

21、*xd0/w(j)*sin(w(j)*t);x=x+xt;end endfor r=1:2subplot(2,1,r) plot(t,x(r,:),grid; xlabel(xxx);ylabel(yyy,num2str(r); end% 把两个分量累加%分别对 x1,x1 绘图1 3。阻尼振动B=0.7;w=5;t=0:0.001:10;x=dsoke(D2x)+2.*B.*(Dx)+w.A2.*x); plot(t,x)7.弹簧振子的物理模型分析本题目用MATLAB处理简单弹簧振子和耦合弹簧振子的物理模型，并对实验数据进行了拟合。(1)简单弹簧振子用劲度系数分别为k1和k2的弹簧，把质量为

22、 m的振子固定在气垫导轨，求振子的速度v,加速度a和位移X。根据牛顿第二定律m*(D2x)=-(k1+k2)*x解微分方程，得运动方程X=A*COS(W*t+ a )对上式微分得速度和加速度v=-A*w*sin(w*t+ a )a=-A*w*w*cos(w*t+ a)其中 w=sqrt(k1+k2)/m);编程syms x m k1 k2;x=dsolve (m*(D2x)+(k1+k2)*x=0,x(0)=0,Dx(0)=10,t);v=diff(x)a=diff(v)(2)耦合弹簧振子用劲度系数为k 1的边弹簧和劲度系数为k2的耦合弹簧连接质量为m1和m2的两振子，求此物理模态的振动系统方

23、程。设x1和x2分别是两振子离开平衡位置的位移，耦合振子的振动方程是 m1*(d(dx1/dt)/dt)=-k1*x1-k2*(x1-x2);m2*(d(dx2/dt)/dt)=-k1*x2-k2*(x2-x1);写成矩阵形式M*(D2x)+K*x=0用eig函数求出矩阵 K- M 的特征值L和特征向量 U, U和L满足U *M*U=I入10 ：L=U *K*U= 02在原始士程 M*x3+K*x=0两端各左乘以 U及在中间的对角矩阵U“U,得U*M* U *U*x +U *K*U *x=0 作变量置换z=z1,z2t= U*x1,x2 t= U*x,z+L*z=0这是一个对角矩阵方程，既可把

24、它分为两个方程z1 +L*z1=0 z2 +2=0 这意味着两种振动模态可以解耦.d* w1= M*与，w2*w2=立*屹，其中w1是第1个模态的固有频率，32是第2个模态的固有频率.由上述的解耦模态中，给出初始条件 x0,xd0,化为z0,zd0求出变量z1,z2再把z转换为x.设速度和 位置的初始条件分别为xd0=xd01,xd02 t, x0=x01,x02 t,则最后得到结果：.x(t)=u1*(u1 t*M*x0*cos( co1*t)+1/w1*u1 t*M*xd0*sin( 1*t)u2 *(u2 t*M*x0*cos( 1*t)+1/w2*u2 t*M*xd0*sin( 2*t

25、)几种不同情况的 MATLAB编程同方向的振动，在t=0时刻，两振子从平衡位置向同一方向移动相同的位移，对此质量相同的两振子来说，两振子以同样的频率，同方向等振幅振动.此时两振子的距离保持常数没有能量的传递，两振子的振动频率仅与两边弹簧有关.即两振子同向振动，此时耦合弹簧不起作用已知条件： 两振子的质量：m1=2;m2=2; 边弹簧的劲度系数：K1=8; 耦合弹簧的劲度系数：K2=4; 两振子的初始位置(在 x轴上)：m1的初始位置：x01=1 ;m2的初始位置：x02=1 两振子的初始速度(以x轴正向取正)：m1的初始速度：x0d1=6;1052 移0爸-5-10-02581时间MATLAB

26、编程clear allm1=.2;m2=.2;K1=8;K2=8;x0=1;0;xd0=6;6;tf=10;M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2;u,L=eig(K,M);t=linspace(0,tf,101);x=zeros(2,101);for s=1:2for j=1:2alfa=sqrt(u(:,s)*M*u(:,s);u(:,s尸 u(:,s)/alfa;m=sqrt(L);m2的初始速度：x0d2=6;%输入各原始参数%初始条件%构成参数矩阵%求广义特征向量和特征值%时间分割和输出变量初始化%分别处理两特征值%解耦后的向量%把特征向量归一化%分别求对

27、应于两特征值的分量w(j)=m(j,j)xt=u(:,j)*(u(:,j)*M*x0*cos(w(j)*t)+u(:,j)*M*xd0/w(j)*sin(w(j)*t);x=x+xt;%把两个分量累加endendfor r=1:2subplot(2,1,r) plot(t,x(r,:),grid;xlabel(时间);ylabel(位移,num2str(r); end运行结果：w =2.00002.8284w =2.00002.8284w =2.00002.8284w =2.00002.8284%分别对x1,x1绘图%绘图振子反方向振动时，在t=0时，使质量相同的两振子位于各自的平衡位置两侧，

28、且到平衡位置的距离相等，此后两振子以相同的频率，方向等幅振动.两振子不发生能量传递.但此时的频率比同方向振动的频率大，周期要短，因为此时耦合弹簧起作用.已知条件：两振子的质量：m1=2;m2=2;边弹簧的劲度系数：K1=8 ;耦合弹簧的劲度系数：K2=4;两振子的初始位置（在x轴上建坐标系）：m1的初始位置：x01=1 ; m2的初始位置：x02=1 ;两振子的初始速度（以 x轴正向取正值）：m1的初始速度：x0d1=6; m2的初始速度：x0d2=-6 ;编程clear all%输入各原始参数%初使条件%构成参数矩阵u（%求广义特征向量和特征值%时间分割和输出变量初始化m1=.2;m2=.2

29、;K1=4;K2=4;x0=1;0;xd0=6;-6;tf=40;M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2;u,L=eig(K,M);t=linspace(0,tf,101);x=zeros(2,101);for s=1:2for j=1:2 alfa=sqrt(u(:,s)*M*u(:,s); u(:,s尸 u(:,s)/alfa;m=sqrt(L);%分别处理两特征值%解耦后的向量%把特征向量归一化%分别求对应于两特征值的分量x=x+xt;endendfor r=1:2subplot(2,1,r) plot(t,x(r,:),grid;xlabel(时间);yla

30、bel(位移,num2str(r); end%把两个分量累加%分别对x1,x1绘图%绘图w(j)=m(j,j) xt=u(:,j)*(u(:,j)*M*x0*cos(wO)*t)+u(:,j)*M*xd0/wO)*sin(wO)*t);运行结果w =2.00002.00002.00002.00002.82842.82842.82842.82847/ f /1/1/ f/ r11TLJr/11/t,i/ 7642-2-4? 0-6 012345678910?.aAf 11/ 1 fT11 J1ri ff,tf/I /642-2-4? 0012345678910?-6 拍振动，在t=0时，两振子到

31、平衡位置的距离不同，振子的质量相同，两振子在耦合弹簧及边弹簧的作用下会发生拍振动.两振子会存在能量传递.当一振子处于平衡位置另一振子离平衡位置距离为A,此两振子在耦合弹簧的作用下会发生拍振动.两振子会存在能量传递.第一个振子振幅变小时，另一个振子振幅增大.已知条件：两振子的质量：m1=2;m2=2;边弹簧的颈度系数：K1=8 ;耦合弹簧的颈度系数:K2=4;两振子的初始位置（在x轴上建坐标系）：m1的初始位置：x01=1 ; m2的初始位置：x02=1 ;两振子的初始速度（以 x 轴正向取正值）：m1的初始速度：x0d1=6; m2的初始速度：x0d2=0 ;编程clear allm1=.2;m2=.2;K1=8;K2=4;x0=1;0;xd0=6;0;tf=40;M=m1,0;0,m2;K=K1+K2,-K2;-K2,K1+K2;%输入各原始参数%初使条件%构成参数矩阵u,L=eig(K,M);t=linspace(0,t