setaspectratio函数_有哪些堪称「数学物理建模软件技巧及神器」，却鲜为人知？…

现在，已经有不少软件可以作物理的动画，进行过程模拟。本教程将结合Mathematica的动画功能，模拟大学物理中两个机械波源所发出波的干涉,以及Matlab在物理学中的应用。

通常，物理教科书上对波的干涉都只讲同频率的两个波源所发出的波的干涉，其实，那只是一个特殊情况。波的干涉是波的叠加性的一种表现，不论两个波源的频率如何，叠加都是存在的，因此，干涉也是存在的。只不过干涉现象的观察手段受到限制，人们只能观察到“稳定的干涉现象”，而同频率的条件可以造成稳定的干涉。如果波源的频率不同，干涉现象就不是稳定的，那么，这时的干涉现象会是怎样的，对此问题的回答可以从干涉的数学表达式上得到，但不够直观。下面，就以在平面上传播的波(如水波)为例，将干涉的数学表达式转化成动画，模拟一般的干涉现象，并观察其特征。

1.波的叠加公式

设在水面上相距为2α的两个地方分别放置半径均为r0的两个柱体，它们各自受到垂直方向的余弦振动源的驱动，周期性地上下起伏，因而引起两个波动向周围传播，并在水面上互相叠加，形成干涉。波源的布置如图1-1所示。

波的振幅为1，各点都作同相位的简谐振动，波的表达式为

其中，

波形示例如下图所示:

下面的程序先画出了两个波源同频率、同相位、同波速振动的波的干涉图，在一个周期内画了21幅干涉图，程序的后面附了其中两幅图，我们可以观察两个时刻干涉图的基本情况，可以看到存在干涉稳定的区域。在程序中，还将产生的每幅图存储在变量animate中，并在最后将其保存为.gif格式的文件，可以用Windows自带IE浏览器等GIF软件播放，以显示动画效果，例如，可以更直观地看到波是如何被波源激发出来，每个区域又如何周期地波动，干涉相消的区域又是如何保持不动的。

将生成的GIF动画保存到d:/animate1. gif,源代码如下:

r0=5.0;a=15.0;v=15.0; ω1=20.0; ω2=20.0;

p={x,y};p1={-a,0};p2={a,0};

animate={ };

φ[t_,x_,y_]:=

Which[ Norm[p-p1] ≤0, Cos[ ω1t],

Norm[p-p2] ≤0, Cos[ω2t], True,

Cos [ω1(t-(Norm[p-p1]-r0)/v)]/Norm[p-p1]+

Cos[ω2(t-(Norm[p-p2]-r0)/v)]/Norm[p-p2]];

Do[

g=Plot3D[φ[t,x,y],{x,-50,50},{y,-50,50},

PlotRange->{-1,1}, PlotPoints->100,

Mesh->False, Boxed->False,

AxesEdge->None, ViewPoint->{2,0,10}];

print[g]; AppendTo [ animate,g],

{t,0,2π/ω1,2π/ω1/20}]

Export[“d:/animate1. gif”, animate]

Clear[“Global’*”]

对以上的模拟说明以下两点。

(1)除了可以画3D图形之外，还可以画波动的密度图，也能展示干涉的情况，但效果稍差一点，我们可以一试。

(2)由于水波离开波源之后振幅与距离的反比关系在减弱，所以，在一个波源的背后，来自另一个波源的波就会比较弱一些，干涉现象就不明显了，基本上是一个波源发出的波在同心地扩散，由此推之，只有在连接两波源的垂线方向(即图中水平方向)，或者两个波源之间的距离与波长相差不多，干涉现象才比较明显。

4.不同波速, 不同频率的干涉

将生成的GIF动画保存到d:/animate3. gif,源代码如下:

r0=2.0;a=15.0;v1=4.0; v2=6.0;ω1=6.0; ω2=10.0;

p={x,y};p1={-a,0};p2={a,0};

animate={ };

φ[t_,x_,y_]:=

Which[ Norm[p-p1] ≤0, Cos[ ω1t],

Norm[p-p2] ≤0, Cos[ω2t], True,

Cos [ω1(t-(Norm[p-p1]-r0)/v1)]/Norm[p-p1]+

Cos[ω2(t-(Norm[p-p2]-r0)/v2)]/Norm[p-p2]];

Do[

g=Plot3D[5*φ[t,x,y],{x,-50,50},{y,-50,50},

PlotRange->{-1,1}, PlotPoints->100,

Mesh->False, Boxed->False,

AxesEdge->None, ViewPoint->{2,0,10}];

print[g]; AppendTo [ animate,g],

{t,0,5×2π/ω1,2π/ω1/10}]

Export[“d:/animate3. gif”, animate]

Print[g], {t,0,5×2π/ω1,2π/ω1/10}]

Clear[“Global’*”]

5.Matlab在电磁学中的应用

电荷在磁场中螺旋线运动的三维动态演示。若电荷以速度v与B成 θ角进入匀强磁场，在这种情况下，可将v分解为与B垂直v1，和与B平行v2的两个分量，

v1=vcosθ，v2=vsinθ，v2不受磁场力的影响，它使电荷沿着B的方向做匀速直线运动；v1会引发磁场力，该磁场力将使电荷在与B垂直的平面内作圆周运动，两种情况的合成使电荷运动的轨迹成为一条螺旋线，如下图所示:

设我舰发射导弹时位置在坐标原点，敌舰在x轴正向d(km)处，其行驶速度为a(km/h)，方向问为x轴夹角为θ,导弹水平飞行线速度为b(km/h)。如上图所示

问题的关键是求出导弹击中敌舰的时间。设t时刻时导弹位置为(x(t),y(t))，那么线速度为

由(1-1)式和(1-2)式得下列微分方程

初始条件x(0)=0,y(0)=0.对于给定的a,b,d, θ,进行计算.当x(t)满足

得,

一般地，应有36

对于测定的d和θ，可用(1-1)式和(1-2)式计算出t.比如d=50，θ=π/2，

用MATLAB写M函数missilefun.m,为了防止分母为0，加了一个小正数1e-8.

missilefun源代码如下:

%M函数missilefun.m,为了防止分母为0，加了1e-8.并且使用附加参数a，b，d，theta传递.

function dy=missilefun(t,y,a,b,d,theta)

dydx=(a*t*sin(theta)-y(2)+1e-8)/…

(abs(d+a*t*cos(theta)-y(1))+1e-8);

dy(1)=b/(1+dydx^2)^0.5;

dy(2)=b/(1+dydx^(-2))^0.5;

dy=dy(:);

计算机模拟动画的方法.对于任意选定的参数a，b，d，θ，T，下面的M函数提供一个导弹追击敌舰演示工具.其中使用了MATLAB动画制作指令getframe和动画播放指令movie.

用MATLAB写M函数missile.m源代码如下:

%M函数missile.m

function m=missile(a,b,d,theta,T)

[t,y]=ode45(@missilefun,[0,T],[0,0],[ ],a,b,d,theta);

x=[d+a*t*cos(theta),a*t*sin(theta)];

n=length(t);

j=n;

for i=1:n

plot(x(i,1),x(i,2),’o’,y(i,1),y(i,2),’r.’);

axis([0 max(x(:,1)) 0 max(x(:,2))]);

hold on;

m(i)=getframe;

if y(i,1)>=x(i,1),

j=i;

break;

end

end

hold off;

movie(m);

legend(‘敌舰’,’导弹’,2);

if j

hold on;

plot(y(j,1),y(j,2),’rh’,’markersize’,18);

hold off;

title([‘导弹将在第’,num2str(t(j)),’小时击中敌舰’]);

else

title([‘导弹在’,num2str(T),’小时内不能击中敌舰’]);

end

对于敌舰速度a=90km/h，导弹速度b=450km/h，距离d=30km，敌舰行驶角度0=0.3m，反应时间T=0.1h，在指令窗口执行

>>missile(90,450,30,0.3*pi,0.1);

操作过程如下:

相关资源：典型自动驾驶仿真软件介绍-iatf16949标准知识培训讲义-交通文档类…