ワンパルスの波を描く

講義の中で衝撃波の話が出てきたのですが,波面の共通接線だけが強めあい,他は打ち消しあう
というのが良くわからなかったので,これはもう計算してしまおうと考えて,以下やってみました.

まずは準備として,ワンパルスの波を作図します.振動数は 1Hz,速度は 4[m/s]です.

以下にプログラムを示します.　周波数は frq = 1 Hzです。波の速度は velo = 4 m/sec です。
波紋の中心は (x0, y0) = ( 0, 0 )です。

プログラムの中で,nami[t,x,y,x0,y0]という関数を定義していますが,これは時刻t=0において
座標(x0,y0)からスタ-トしたワンパルスの波が,時刻tのとき座標(x,y)においてどのような変位で
あるかを与えるものです.

nami[t_,x_,y_,x0_,y0_]:=
           If[(velo*(t-shuki)<r[x,y,x0,y0]) && (r[x,y,x0,y0]<velo*t),
              Sin[-2*Pi*(r[x,y,x0,y0]-velo*(t-shuki))/hachou]*warur[x,y,x0,y0],0];

ワンパルスにするためにIf文を用いて波を切り抜いています.

波の開始地点からの距離を求めるのに r[x,y,x0,y0]という関数も定義しています.

r[x_,y_,x0_,y0_]:=Sqrt[(x-x0)^2 + (y-y0)^2];

波の振幅は伝わるにつれて,1/rに比例して小さくなっていきます.そこで 1/rを掛け算することが
必要になりますが,この値が無限大になってしまうのが怖いので、rで割り算する関数 warur[x,y,x0,y0]も
用意しました.r=0のところで多少細工をしています.

warur[x_,y_,x0_,y0_]:=If[Sqrt[(x-x0)^2 + (y-y0)^2]==0,1,
　　　　　　　　　　　　(1-Exp[-Sqrt[(x-x0)^2 + (y-y0)^2]])/Sqrt[(x-x0)^2 + (y-y0)^2]];

作図は3秒経過した後の様子を作図しました.以下プログラムです。

Remove["Global `*"];
Timing[
  frq=1;
  shuki=1/frq;
  r[x_,y_,x0_,y0_]:=Sqrt[(x-x0)^2 + (y-y0)^2];
warur[x_,y_,x0_,y0_]:=If[Sqrt[(x-x0)^2 + (y-y0)^2]==0,1,
  　　　　　　　　　(1-Exp[-Sqrt[(x-x0)^2 + (y-y0)^2]])/Sqrt[(x-x0)^2 + (y-y0)^2]];
  
  velo=4;
  hachou=velo/frq;
  
nami[t_,x_,y_,x0_,y0_]:=If[(velo*(t-shuki)<r[x,y,x0,y0]) && (r[x,y,x0,y0]<velo*t),
  　　　　　　Sin[-2*Pi*(r[x,y,x0,y0]-velo*(t-shuki))/hachou]*warur[x,y,x0,y0],0];
  
  t=3;
  x0=0;
  y0=0;
  
  Plot3D[nami[t,x,y,x0,y0],{x,-20,20},{y,-20,20},PlotRange → {-0.5,0.5},PlotPoints →80];
  ]

せっかくですから,この絵を動かして見ましょう.時刻を変えながらたくさんの絵を描くだけです.
    Do文を用いて繰り返しています.51枚のグラフを描くのですが,Core2Duoの2.4GHzで8秒ぐらいかかります.
    メモリ-も必要です.以下プログラムです.

Remove["Global `*"];
Timing[
  frq=1;
  shuki=1/frq;
  r[x_,y_,x0_,y0_]:=Sqrt[(x-x0)^2 + (y-y0)^2];
  warur[x_,y_,x0_,y0_]:=If[Sqrt[(x-x0)^2 + (y-y0)^2]==0,1,
  　　　　　　　　　　　　　(1-Exp[-Sqrt[(x-x0)^2 + (y-y0)^2]])/Sqrt[(x-x0)^2 + (y-y0)^2]];
  
  velo=4;
  hachou=velo/frq;
  
  nami[t_,x_,y_,x0_,y0_]:=If[(velo*(t-shuki)<r[x,y,x0,y0]) && (r[x,y,x0,y0]<velo*t),Sin[-2*Pi*(r[x,y,x0,y0]-velo*(t-shuki))/hachou]*warur[x,y,x0,y0],0];
  
  x0=0;
  y0=0;
  
  Do[
    t=(8/50)*i;
    Plot3D[nami[t,x,y,x0,y0],{x,-20,20},{y,-20,20},PlotRange → {-0.5,0.5},PlotPoints →80];
    ,{i,0,50,1}];
  ]

グラフを描いたら,そのセルをすべて選択して,メニュ-にある[セル]-[グラフィックスのアニメ-ション化]で
動かすことができます(動画をお見せしたいところですが,少しお待ちください).

この動画をaviファイルとして保存します.これには Exportという命令を使います.

Export["test.avi",gsnaplist,"AVI"];

動画の入れ物として gsnaplistというリストを用意しています.

gsnaplist=Table[Null,{i,60+2}];

このリストに描いたグラフを収納して,
後は Exportを用いるだけです.

Remove["Global `*"];
Timing[
  frq=1;
  shuki=1/frq;
  r[x_,y_,x0_,y0_]:=Sqrt[(x-x0)^2 + (y-y0)^2];
  warur[x_,y_,x0_,y0_]:=If[Sqrt[(x-x0)^2 + (y-y0)^2]==0,1,
  　　　　　　　　　　　　　　　(1-Exp[-Sqrt[(x-x0)^2 + (y-y0)^2]])/Sqrt[(x-x0)^2 + (y-y0)^2]];
  
  velo=4;
  hachou=velo/frq;
  
  nami[t_,x_,y_,x0_,y0_]:=If[(velo*(t-shuki)<r[x,y,x0,y0]) && (r[x,y,x0,y0]<velo*t),Sin[-2*Pi*(r[x,y,x0,y0]-velo*(t-shuki))/hachou]*warur[x,y,x0,y0],0];
  
  x0=0;
  y0=0;
  
gsnaplist=Table[Null,{i,60+2}];
  
  Do[
    t=(8/60)*i;
    gsnaplist[[i+1]]=Plot3D[nami[t,x,y,x0,y0],{x,-20,20},{y,-20,20},
    　　　　　　　　　　PlotRange → {-0.5,0.5},PlotPoints →80,DisplayFunction->Identity ];
    ,{i,0,60,1}];
  
  SetDirectory["d:\\today"];
  Export["test.avi",gsnaplist,"AVI"];
  ]

これで作成されるaviファイルは無圧縮のものです.圧縮にはwindows media encoder が良いと思いますが,
aviファイルの縦横のサイズによって上手くいかないことがあります.

環境設定のグラフィックの設定で,グラフィックの境界ボックスのうち,ImageSizeを720x720 にすると
上手く変換できました(理由は良くわかりません).

このようにしてこの動画 ファイルができました.

このページで使用したプログラム


　プログラム１ 　プログラム２ 　プログラム３

Created by Mathematica  (May 31, 2007)