物理と数学シリ-ズ　With Mathematica

§３　今から未来へ　（微分方程式基礎）

本講座の受講の前に次の講座を修了していることが理想である。
Mathematicaの基礎
１．物理と数学シリーズ　微分と積分
２．物理と数学シリーズ  三角関数の基礎

前回微分について学んだ微分はある関数の傾きを表している。これはいいかえると次の状態に移る時、どれくらい変化するかを示している。例えば y=2x という関数を考えるとこの傾きは２だからｘが１増えればｙは２増える。従って次のような簡単な関係がｙの値と微分した値にはある。

yの変化量 =[微分係数] × xの変化量 (1)

普通の関数ではｘの変化量の値はどんな数でも選べられるし、ある小さな数を選んだとしてもその数よりも小さい数を選ぶことが必ずできる。そこでｘの変化量はある一定な値であると決めよう。そうすればｘがどれだけ変化したかはこの基本的な値の整数倍に必ずなる。
例えば先のｙ＝２ｘではx=1の時yは２、ｘ=2の時はｙは４になりｘが１だけ変化すればちょうどｙは微分係数だけ変化する。

これから次のように書き直すことができる。

yのn番目の値 =[定数] × yの (n - 1) 番目の値 (2)

では具体的に例を見ていこう。

例えば人口がどれだけ増えていくかを考える。未来の人口は現在の人口だけで決まるとしよう。人口をｙとして、微分係数を定数Rで表すと次のようになる。

Out[74]=

y[n] = R y[n - 1] (3)

【Mathematica】

In[38]:=

y[n_] := R y[n - 1]

y[0] = 1 ;

R = 1.1 ;

tb = Table[y[n], {n, 1, 50}] ;

ListPlot[tb]

[Graphics:HTMLFiles/index_9.gif]

Rが1.1ということは今より次の世代はわずか１０％増えることになる。しかし、上のグラフを見てわかるように時がたつと急激に人口は増えていく。
次に同じ関係を次の式で表す。
y'でyのｘについての微分を表すとすると式(3)に変わる式は次のような微分を含んだ式になる。

Out[74]=

y = R y ' (4)

微分を含まない式と微分を含んだ式が混在して一つの方程式になっているがこれを微分方程式という。
これをMathematicaのDSoleveを用いて解いてみよう。yをｘの関数と明示するためにf(x)で表すことにしよう。

In[5]:=

R = 1.1 ;

g1 = f[x]/.DSolve[{f[x] == R f '[x], f[0] == 1}, f[x], x] ;

fg = Part[g, 1]

Plot[fg, {x, 0, 5}]

Out[7]=

E^(0.909091 x)

[Graphics:HTMLFiles/index_16.gif]

期待したとおり指数関数が出てきてグラフは先の点のグラフとほとんど一致している。
微分方程式を解く場合は初期値が重要でこの場合はｘが０の時１として解いている。簡単のための数字から出発してるので一人からどうやって人口を増やすかと悩まないように！

前回の微分の基礎で見たように高校での運動方程式は微分方程式である。例えば次ぎのように摩擦の無い水平面で質量2kgの物体を水平に10Nの力Fで引っ張る場合の運動方程式は次のようになった。

[Graphics:HTMLFiles/index_17.gif]

Out[74]=

m a = F (5)

これは微分方程式では次のようになる。だたし、今度はx'はｘのｔについての微分を表す。従ってx'は速さｖ、x''は加速度aを表す。

Out[74]=

m x'' = F (6)
x'' = F/m (6)

【Q1】　Fを10、ｍを２としてこの式をMathematicaで解いてみよ。

【Mathematica】

In[23]:=

Clear[x]

g1 = x[t]/.DSolve[{x''[t] == 25, x[0] == 0, x '[0] == 0}, x[t], t] ;

fg = Part[g, 1]

Plot[fg, {t, 0, 1}]

Out[25]=

(25 t^2)/2

[Graphics:HTMLFiles/index_26.gif]

上の問題の初期条件は時刻０で原点の位置x[0]=0、また時刻０で初速０ x'[0]=0としている。
結果はおなじみの公式が出てきてｔの2次関数となる。

【Q2】　先の問題を動摩擦係数0.5としてMathematicaで解いてみよ。

In[32]:=

Clear[x] ; Clear[g] ;

g = 9.8 ;

g1 = x[t]/.DSolve[{x''[t] == 25 - 0.5 g, x[0] == 0, x '[0] == 0}, x[t], t] ;

fg = Part[g1, 1]

Plot[fg, {t, 0, 1}]

Out[35]=

10.05 t^2

[Graphics:HTMLFiles/index_33.gif]

Out[36]=

- Graphics -

摩擦が入ったとしても微分方程式の形としては基本は次のように
x''[t]=定数
変わっていない。そのため結果はｔの2次関数で傾きは減ったがグラフの形はおよそ同じになる。

それでは次に自由落下を考えてみよう。ただし、空気抵抗として速度に比例した抵抗を考える。つまり物体の落下速度が速くなればなるほどそれに比例して空気抵抗が働くとする。この時の運動方程式は次のようになる。

Out[74]=

m x'' = m g - k v (7)

【Q3】　上の運動方程式質量1kg、ｋを0．5としてMathematicaで解いてみよ。

【Mathematica】

In[67]:=

Clear[x] ; Clear[g] ; Clear[g1] ; Clear[y] ;

g = 9.8 ;

g1 = y[t]/.DSolve[{y''[t] == g - 0.5y '[t], y[0] == 0, y '[0] == 0}, y[t], t] ;

fg = Part[g1, 1]

Plot[fg, {t, 0, 1}]

Out[70]=

2.71828^(-0.5 t) (39.2 - 39.2 2.71828^(0.5 t) + 19.6 2.71828^(0.5 t) t)

[Graphics:HTMLFiles/index_42.gif]

Out[71]=

- Graphics -

なぁんださっきのグラフとかわらないじゃん！と言ってはいけない。さっきのグラフは等加速運動なわけだから速さは永遠に増えていく。
しかし、このグラフは縦軸をｙ（落下距離）としていたらその様子は見えにくい。では縦軸は速さ、すなわちy'[t]のグラフを描いてみればいい。
Mathematicaでは次のようになる。

In[90]:=

Clear[x] ; Clear[g] ; Clear[g1] ; Clear[y] ;

g = 9.8 ;

g1 = y[t]/.DSolve[{y''[t] == g - 0.5y '[t], y[0] == 0, y '[0] == 0}, y[t], t] ;

fg = Part[g1, 1]

fv = D[fg, t]

Plot[fv, {t, 0, 7}]

Out[93]=

2.71828^(-0.5 t) (39.2 - 39.2 2.71828^(0.5 t) + 19.6 2.71828^(0.5 t) t)

Out[94]=

2.71828^(-0.5 t) (0. 2.71828^(0.5 t) + 9.8 2.71828^(0.5 t) t) - 0.5 2.71828^(-0.5 t) (39.2 - 39.2 2.71828^(0.5 t) + 19.6 2.71828^(0.5 t) t)

[Graphics:HTMLFiles/index_52.gif]

Out[95]=

- Graphics -

グラフは時間がたつと速さが一定になりそれ以上増えないことを表す。実際にある物体が自由落下するときある一定な速さに近づいていく。

これは最初の人口の問題にもどると増加の割合が負になる場合であるから次のように考えることができる。

In[561]:=

Clear[y] ;

y[n_] := g + k y[n - 1]

g = 1 ;

y[0] = 0 ;

k = 0.9 ;

tb = Table[y[n], {n, 0, 40}] ;

ListPlot[tb, PlotJoined -> True]

[Graphics:HTMLFiles/index_61.gif]

Out[567]=

- Graphics -

■単振動

質量ｍのおもりにバネ定数ｋのバネをつけて振動させる場合に単振動の一般式は次のように表すことができた。

m a == -k x

加速度aは位置ｘの2回微分として書きなおすと

x''[t] + ω^2 x[t] == 0

ただし、

ω^2 = k/m

である。これに初期条件としてｔ＝０でｖ＝０，ｘ＝Aとすれば
次のようにこの解を求めさせることができる。

In[1]:=

DSolve[{x''[t] + ω^2 x[t] == 0, x '[0] == 0, x[0] == A}, x[t], t] ;

y = x[t] /. Part[%, 1]

Plot[y /.{A -> 1, ω -> 1}, {t, 0, 4π}]

Out[2]=

A Cos[t ω]

[Graphics:HTMLFiles/index_70.gif]

Out[3]=

- Graphics -

Out[24]=

A Cos[t ω] (8)

単振動は円運動の射影であった。従って射影の際に情報が落ちているので単振動１つから逆に円運動に戻すことはできない。
しかし、次のように2次元の座標を用意し、ｘ座標に式（８）の結果を入れ、ｙ座標にはこれと π/2 だけずらした値をいれｔをちょうど１周期だけ変化させてみる。
ただし、簡単のため ω=1,A=1 とする。

In[37]:=

ParametricPlot[{Cos[t], Sin[t]}, {t, 0, 2π}, AspectRatio -> Automatic] ;

[Graphics:HTMLFiles/index_74.gif]

円が再現された。
さらに少し進んで ωの値をｘ座標、Ｙ座標で変化させたらどうなるか
次のようにして試してみよう。

In[10]:=

ParametricPlot[{Cos[5t], Sin[3t]}, {t, 0, 2π}, AspectRatio -> Automatic]

[Graphics:HTMLFiles/index_76.gif]

円ではないが連続した周期的な図形が描かれた。

■空気抵抗

空気抵抗が無い場合の運動方程式は次のようになった。

m x''[t] = m g

これを解くには次のようにする。自由落下としてｔ＝０でｘ＝０、ｖ＝０とすると

In[84]:=

Clear[g] ; Clear[y1] ;

DSolve[{x''[t] == g, x '[0] == 0, x[0] == 0}, x[t], t] ;

y1 = x[t] /. Part[%, 1]

Plot[y1 /.g -> 9.8, {t, 0, 2}, AspectRatio -> 1]

Out[86]=

(g t^2)/2

[Graphics:HTMLFiles/index_83.gif]

Out[87]=

- Graphics -

予想通りこれは放物線であり等加速度運動である。

次に空気抵抗がありこれが速度に比例する抵抗を持つ場合の運動方程式は比例定数をaとして

m x''[t] = m g - a x '[t]

ここで新に比例定数ｋを

k = a/m

とおけば次のようにして解くことができる。

In[88]:=

Clear[g] ; Clear[y2] ;

DSolve[{x''[t] == g - k x '[t], x '[0] == 0, x[0] == 0}, x[t], t] ;

y2 = x[t] /. Part[%, 1]

Plot[y2 /.{g-> 9.8, k -> 1}, {t, 0, 1}, AspectRatio -> 1]

Out[90]=

E^(-t) g (1 - E^t + E^t t)

[Graphics:HTMLFiles/index_92.gif]

Out[91]=

- Graphics -

これを微分して速さのグラフV-tグラフを描かせてみよう。自由落下の場合のy1を微分したときと比べてみよう。

In[92]:=

g1 = D[y1, t] /.{g -> 9.8, k -> 1} ;

g2 = D[y2, t] /.{g -> 9.8, k -> 1} ;

Plot[{g1, g2} , {t, 0, 1}]

[Graphics:HTMLFiles/index_97.gif]

Out[94]=

- Graphics -

tが小さい時２つのグラフは一致しているがｔが大きくなると速さに差がどんどんできてくるのがわかる。

この場合は重力が働いているのでどちらにしろ物体はどこまでも落ちていくし、速度は抵抗があるからといって０になるわけではなくある値に近づいていく。
では速度に比例した抵抗が水平面に働くとして水平面上に質量ｍの物体を初速V0で滑らす場合はどうか？この時の運動方程式は次のようになる。

m x''[t] = -a x '[t]

ここで新に比例定数hを先と同じように

h = a/m

とおけば次のようにして解くことができる。
ただし、初期条件として
ｔ＝０でｖ＝V0、ｘ＝０とする。

In[95]:=

DSolve[{x''[t] == -h x '[t], x '[0] == V0, x[0] == 0}, x[t], t] ;

y3 = x[t] /. Part[%, 1]

Plot[y3/.{g-> 9.8, h -> 0.5, V0 -> 1}, {t, 0, 4}, AspectRatio -> 1]

Out[96]=

(E^(-h t) (-1 + E^(h t)) V0)/h

[Graphics:HTMLFiles/index_105.gif]

Out[97]=

- Graphics -

これを微分して速さのグラフV-tグラフを描かせてみよう。

In[98]:=

g3 = D[y3, t] /.{g -> 9.8, h -> 0.5, V0 -> 1}

Plot[g3, {t, 0, 4}, AspectRatio -> 1]

Out[98]=

1 - E^(-0.5 t) (-1 + E^(0.5 t))

[Graphics:HTMLFiles/index_110.gif]

Out[99]=

- Graphics -

予想通り時間がたつほど速さは０に近づく。
ここで再び最初の単振動のV-tグラフを描かせてみよう。

In[100]:=

g = D[y, t] /.{ω -> 1, A -> 1}

Plot[g, {t, 0, 4π}, AspectRatio -> 1]

Out[100]=

-Sin[t]

[Graphics:HTMLFiles/index_115.gif]

Out[101]=

- Graphics -

速度が０からマイナス方向に増えていってるのが理解できる。これはｘーｔグラフと対比させると物体の位置と速度の関係が微分で結ばれさらに位相差がπ/2 だけあることも次のようにグラフを重ねれば理解できる。

In[102]:=

y = y/.{ω -> 1, A -> 1}

Plot[{y, g}, {t, 0, 4π}, AspectRatio -> 1]

Out[102]=

Cos[t]

[Graphics:HTMLFiles/index_121.gif]

Out[103]=

- Graphics -

ではこの単振動に先の空気抵抗が入った場合はどうなるか考えてみよう。
運動方程式は簡単に次のようになると予想される。

m x''[t] = -a x '[t] - k x[t]

ここで新に比例定数h,ωを先と同じように

ω^2 = k/m, h = a/m

とおけば

x''[t] + h x '[t] + ω^2 x[t] == 0

とおけば次のようにして解くことができる。
ただし、初期条件として
ｔ＝０でｖ＝0、x=Aとする。

In[104]:=

DSolve[{x''[t] + h x '[t] + ω^2 x[t] == 0, x '[0] == 0, x[0] == A}, x[t], t] ;

yy = x[t] /. Part[%, 1]

Plot[yy/.{A -> 1, ω -> 1, h -> 0.5}, {t, 0, 4π}]

Out[105]=

1/(2 (h^2 - 4 ω^2)^(1/2)) (A (-E^(1/2 t (-h - (h^2 - 4 ω^2)^(1/2))) h + E^(1/2 t (-h + (h^2 ... 4 ω^2)^(1/2))) (h^2 - 4 ω^2)^(1/2) + E^(1/2 t (-h + (h^2 - 4 ω^2)^(1/2))) (h^2 - 4 ω^2)^(1/2)))

[Graphics:HTMLFiles/index_130.gif]

Out[106]=

- Graphics -

減衰してやがて静止するであろうことがわかる。
ではこの場合の速度はどうなるか微分して見てみよう。

In[107]:=

gg = D[yy, t]/.{A -> 1, ω -> 1, h -> 0.5} ;

Plot[gg/.{A -> 1, ω -> 1, h -> 0.5}, {t, 0, 4π}]

[Graphics:HTMLFiles/index_134.gif]

Out[108]=

- Graphics -

位置とはやはり π/2だけずれて減衰していく周期は抵抗がない時と変化していないこともわかる。
g3とyyのグラフと共に表してみよう。

In[109]:=

yy = yy/.{A -> 1, ω -> 1, h -> 0.5}

Plot[{g3, yy}, {t, 0, 4π}]

Out[109]=

(0. - 0.258199 I) ((-0.5 + 1.93649 I) E^((-0.25 - 0.968246 I) t) + (0.5 + 1.93649 I) E^((-0.25 + 0.968246 I) t))

[Graphics:HTMLFiles/index_140.gif]

Out[110]=

- Graphics -

■　アニメーション
次に下図のような単振動の様子をMathematicaを用いて実際にその運動を作って動かしてみよう。
Ｑ　まずｔ＝０で原点を通り振幅がＡのバネの運動方程式を作成し、解いてみよ。

In[245]:=

DSolve[{x''[t] + ω^2 x[t] == 0, x '[0] == ω A, x[0] == 0}, x[t], t] ;

g1 = x[t] /. Part[%, 1]

Plot[g1 /.{A -> 1, ω -> 1}, {t, 0, 4π}]

Out[246]=

Sin[t]

[Graphics:HTMLFiles/index_146.gif]

Out[247]=

- Graphics -

上記の例では原点での速さが Aωであることを利用している。

Out[118]=

A Sin[t ω]

[Graphics:HTMLFiles/index_149.gif]

Out[119]=

- Graphics -

In[248]:=

A=1;
ω=1;
Plot[g1,{t,-2π,2π}]

[Graphics:HTMLFiles/index_151.gif]

Out[250]=

- Graphics -

動画のプログラムは次のようになる。Ｐａｒｔコマンドを利用して、
実際に絵を動かすときは下の例では２４枚の絵のセルを全て選択してCtrL＋Yキーを押す。

In[252]:=

rg = 2π ;

n = 24 ;

h = 4/5 ;

Tb = Table[g1, {t, -rg, rg, rg/n}] ;

ra = Max[Tb] ;

rb = 0.1 ;

gh = {RGBColor[h, h, h], Rectangle[{0, 0}, {2 (ra + rb), 2 (ra + rb)}]} ;

Do[ gp = {RGBColor[0, 0, 1], Disk[{(ra + rb), Part[Tb, i] + (ra + rb)}, rb]} ; Show[Graphics[{gh, gp}], AspectRatio -> 1], {i, 1, n}]

[Graphics:HTMLFiles/index_185.gif]

減衰振動のアニメーションも同様に次ぎのようにして作成することができる。

In[260]:=

a=2;
b=1/2;
DSolve[{y2''[x]+a y2[x]+b y2'[x]==0,y2'[0]==10,y2[0]==0},y2[x],x]
g2=y2[x] /.%
Plot[g2,{x,-5,5}]

Out[262]=

{{(2.71828^(-0.9 t) (1.23457 - 1.23457 2.71828^(0.9 t) + 1.11111 2.71828^(0.9 t) t))[x] -> (40 E^(-x/4) Sin[(31^(1/2) x)/4])/31^(1/2)}}

Out[263]=

{(40 E^(-x/4) Sin[(31^(1/2) x)/4])/31^(1/2)}

[Graphics:HTMLFiles/index_188.gif]

Out[264]=

- Graphics -

In[305]:=

rg = 4π ;

n = 24 ;

h = 4/5 ;

Tb = Table[Part[g2, 1], {x, 0, rg, rg/n}] ;

ra = Max[Tb] ;

rb = 0.3 ;

gh = {RGBColor[h, h, h], Rectangle[{0, 0}, {2 (ra + rb), 2 (ra + rb)}]} ;

Do[ gp = {RGBColor[1, 0, 0], Disk[{(ra + rb), Part[Tb, i] + (ra + rb)}, rb]} ; Show[Graphics[{gh, gp}], AspectRatio -> 1], {i, 1, n}]

[Graphics:HTMLFiles/index_222.gif]

Created by Mathematica  (July 8, 2007) Valid XHTML 1.1!