23章アニメーション

これまでの学習をもとに、さまざまなアニメーションを作成してみましょう。

1. 波の伝播

>with(plots):

>animate(sin(x-t),x=0..4*Pi,t=0..10)

と入力し、描かれたグラフをクリックして表示されるメニューバー上のアニメーション再生ボタンをクリックすると、アニメーションが始まります。メニューバー上のボタンを使って、アニメーションの実行速度を変えたり連続再生モードに設定することなどができます。

注意 u(x,t)=sin(x-t)は微分方程式u_tt=u_xx(波動方程式)を満たします。すなわち波動方程式の解です。

　進行波解とよばれる微分方程式の解のシュミレーションを行ったり、らせんを一定の角速度で回転させてみましょう。

>animate(tanh(x-t),x=-20..20,t=0..25)
(進行波解のシミュレーション)

>animate([r*cos(r-t),r*sin(r-t),r=0..6*Pi],t=0..10,numpoints=100,scaling=constrained)
（回転するらせんのシミュレーション）

注意
(1)u(x,t)=tanh(x-t)は微分方程式u_t　u_xx+(2u-1)(1-u²)を満たします。このように、ある関数φ(y)と定数cを用いてu(x,t)=φ(x-ct)と表せる一定波形で一定速度cで進む解を進行波解といいます。数理生物学や燃焼理論への応用上も重要な解です。
(2)曲線(rcosr,rsinr)、r∈[a,b]はアルキメデスのらせんとよばれます。例題の回転するらせんは、結晶成長を記述する数理モデルに現れる解です。

Kdv(Korteweg-de Iries) 方程式　u_t+6uu_x+u_xxx=0のソリトン解(孤立波解) u(x,t)=2k²sech²(kx-4k³t)
(kは定数)や２ソリトン解 u(x,t)=2(log(1+ce^2kx-8k³t+de^2hx-8h³t+cd(k-h)²e^{2kx-8k³t+2hx-8h³t}/(k+h)²))_xx
(c,d,k,hは定数)のシミュレーションを行ってみましょう。

>animate(2*sech(x-4*t)²,x=-25..25,t=-5..5 ,numpoints=200)
（ソリトン解(k=1の場合）のシミュレーション)

>v:=(x,t)->1+exp(2*x-8*t)+exp(x-t)+9*exp(2*x-8*t+2*x-t)/16

>u:=(x,t)->2*diff(log(v(x,t)),x,x)

>animate(u(x,t),x=-30..30,t=-6..6,numpoints=200)
(2ソリトン解(c=d=k=1,h=1/2の場合)のシミュレーション）

2. サイクロイド
サイクロイドは、円が x 軸に接しながら回転するとき、その周上に固定された点 P の軌跡として定義されます。円の半径と a とすると、

x = a (θ-sin θ)
y = a (1-cos θ)
と表せます。a=1 の場合に、θを 0 から 2π まで動かした時の軌跡を描いてみましょう。

>with(plots):

>animatecurve([t-sin(t), 1-cos(t), t = 0 .. 2*Pi], scaling = constrained, frames = 16)
(サイクロイドを描く）

>animate([cos(s)+t, sin(s)+1, s = 0 .. 2*Pi], t = 0 .. 2*Pi, scaling = constrained)
(回転する円のシミュレーション）

>animate([s*(t-sin(t))+(1-s)*t, s*(1-cos(t))+1-s, s = 0 .. 1], t = 0 .. 2*Pi, scaling = constrained)
(動径のシミュレーション）

>display({%, %%, %%%})

3. 振り子の運動
原点 (0,0) に回転軸をもつ振り子が、重力の作用により運動する様子を考えます。
時刻 t において、振り子の糸が鉛直下向きと半時計回りを正の向きにとって、s(t) の角度をなす時

ml s''(t) = -mg sin s(t)
が成り立ちます。ここで、m はおもりの質量、g は重力加速度、l は糸の長さを表します。 l=g の場合に、初期条件 s(0)=0、s'(0)=0.8 を満たす解を求め、シミュレーションしてみましょう。厳密解は求められないので、dsolve コマンドのオプションで numeric と指定して、数値解を求めることにします。

>with(plots):

>dsolve({diff(s(t), t, t) = -sin(s(t)), s(0) = 0, (D(s))(0) = .8}, numeric, output = listprocedure);

>odeplot(%, [sin(s(t)), -cos(s(t))], t = 0 .. 30, style = point, frames = 100, scaling = constrained)

実習23.1　上の振り子の運動について、初期条件が以下の場合についてのアニメーションを作成してみましょう。
　　(1) s(0)=π/2、s'(0)=0
　　(2) s(0)=0、s'(0)=1.8
　　(3) s(0)=0、s'(0)=2.2

[正解例]

23章 アニメーション

23章アニメーション