23章 アニメーション


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

1. 波の伝播

>with(plots):

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


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

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

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

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

Kdv(Korteweg-de Iries) 方程式　ut+6uux+uxxx=0のソリトン解(孤立波解)
u(x,t)=2k2sech2(kx-4k3t)

(kは定数)や２ソリトン解
u(x,t)=2(log(1+ce2kx-8k3t+de2hx-8h3t+cd(k-h)2e2kx-8k3t+2hx-8h3t/(k+h)2))xx

(c,d,k,hは定数)のシミュレーションを行ってみましょう。

>animate(2*sech(x-4*t)2,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

[正解例]