宇宙の歴史を俯瞰することは人間にはできません。我々にできるのは、今、この一瞬に何が起きているかを観察することです。
ハレー彗星の動きを、太陽を一周するまで観測することは不可能です(そんな長生きの人はいないでしょう)。
我々が観測できるのは、多くて数年です。
宇宙から見ると一瞬にすぎない我々の観察から、宇宙全体の謎を解き明かす方法、それが微分方程式です。
Ⅰ 宇宙船の運行
速度や加速度は、瞬間を観測することにより得られます。速度や加速度の関係を式で表したものが、運動の方程式です。
速度は座標の1階微分、加速度は2階微分ですから、運動の方程式は微分方程式にほかなりません。
身近な運動の方程式を、微分方程式として捉え、これを解くことで質点の運動を解析しましょう。
高校の範囲で解ける微分方程式は限られています。ここでは、数値計算法を用いて、運動の謎に迫ります。
あなたは宇宙船に閉じ込められています。あなたが見えるのは、宇宙船の速度計のみです。手元に地図があり、進んでいる道はわかるのですが、外が見えないので、どこまで来ているのかわかりません。
①次のデータは、速度の情報から位置を算出したものです。(赤の計算結果を確かめましょう)
時刻(秒) | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
速度(m) | 0 | 2 | 4 | 6 | 8 | 10 | 12 | 14 | 16 | 18 | 20 | 22 | 24 |
位置計算 | 0 | 2 | 6 | 12 | 20 | 30 | 42 | 56 | 72 | 90 | 110 | 132 | 156 |
② 算出した位置をプロットします。
これは、どんな関数のグラフに近いでしょうか?
赤のグラフは理論上の数値です。
速度 \(v(t)=2t\) であるので、理論上の位置は \(x(t)=t^2\) となります。
このように近似を用いた数値計算では誤差が生じます。
Ⅱ バネと単振動
バネの伸びと力のあいだには、フックの法則が成立します。
バネが引く力は、バネの伸びに比例するというものです。
式で表すと、
\(f=-kx\) ( \(k\) はバネ定数)
運動の方程式は \(f=ma\)( \(a\) は加速度、 \(m\) はおもりの質量)であるので
\(ma=-kx\)
簡単のため、初期条件を \(m=1\:\:\rm[kg]\:,\:k=1\:\:\rm[N/m]\) とすると
\(a=-x\)
加速度 \(a\) 、おもりの位置 \(x\) は時刻 \(t\) の関数なので、\(a=a(t)\:,\:x=x(t)\) とすると
\(a(t)=-x(t)\)
加速度 \(a\) は位置 \(x\) の2階微分であるので \(a(t)=\dfrac{d^2x}{dt^2}\)
よって、微分方程式 \(\dfrac{d^2x}{dt^2}=-x(t)\) が成立します。
さて、バネを \(0.1\rm \:[m]\) 伸ばして手を放すとします。
\(x(0)=0.1 \rm \:\:[m]\) 、速度:\(v(0)=0\:\:\rm[m/s]\)
①速度、位置、加速度の順に各段に記入し位置計算をしました。
\(a(t)=-x(t)\) の式を用いると、
\(x(0)=0.1\) なので \(a(0)=-0.1\)
0.5秒間隔なので、
\(v(0.5)=0+0.5\times (-0.1)=-0.5\)
\(x(0.5)=0.1+0.5\times (-0.5)=0.75\) \(\therefore a(0.5)=-0.75\)
(電卓を用いて、赤の計算結果を確かめましょう)
時刻 t |
位置 x(t) |
速度 v(t) |
加速度 a(t) |
0.0 | 0.100 | 0.000 | -0.100 |
0.5 | 0.075 | -0.050 | -0.075 |
1.0 | 0.031 | -0.088 | -0.031 |
1.5 | -0.020 | -0.103 | 0.020 |
2.0 | -0.067 | -0.093 | 0.067 |
2.5 | -0.097 | -0.060 | 0.097 |
3.0 | -0.102 | -0.011 | 0.102 |
3.5 | -0.082 | 0.040 | 0.082 |
4.0 | -0.042 | 0.081 | 0.042 |
4.5 | 0.009 | 0.102 | -0.009 |
5.0 | 0.058 | 0.097 | -0.058 |
5.5 | 0.092 | 0.068 | -0.092 |
6.0 | 0.103 | 0.022 | -0.103 |
6.5 | 0.089 | -0.029 | -0.089 |
7.0 | 0.052 | -0.074 | -0.052 |
7.5 | 0.002 | -0.099 | -0.002 |
8.0 | -0.048 | -0.101 | 0.048 |
②計算した位置をプロットすると、
これは、どんな関数のグラフに近いでしょうか?
バネの運動は単振動と呼ばれ三角関数で表されます。
Ⅲ コーヒーは冷めないうちに
コーヒーの冷める速さは、そのときのコーヒーと室温の温度差に比例します。
室温を20℃、 \(t\) 分後のコーヒーの温度を \(x(t)\) とするとき、
\(\dfrac{dx}{dt}=-\alpha(x-20)\)
比例定数 \(\alpha=0.08\) 、初めのコーヒーの温度を80℃として、10分後の温度を算出しましょう。
①コーヒーの温度、室温との差の順に各段に記入しましょう。
\(x(1)=80-0.08\times 60=75.2\) 室温との差は \(75.2-20=55.2\)
(赤の計算結果を確かめましょう)
時刻(分) | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | \(\cdots\) |
お湯の温度 | 80.0 | 75.2 | 70.8 | 66.7 | 63.0 | 59.5 | 56.4 | 53.5 | 50.8 | 48.3 | 46.1 | \(\cdots\) |
室温との差 | 60.0 | 55.2 | 50.8 | 46.7 | 43.0 | 39.5 | 36.4 | 33.5 | 30.8 | 28.3 | 26.1 | \(\cdots\) |
②計算した温度をプロットします。
これは、どんな関数のグラフに近いでしょうか?
理論上は \(x=60e^{-0.08t}+20\) ( \(t\) の単位は分のため、60倍がつく) となります。
Ⅳ マルサスの人口論
マルサス(1766~1834)はイギリスの経済学者です。「出生数と死亡数は人口の大きさに比例する」
「人口の増加は人口の大きさに比例する」という仮説のもとで、人口論を展開しました。
日本の人口の変化を数学的に予測しましょう。人口の増加率は現在の人口に比例します。
時刻 \(t\) における人口を \(x(t)\) とすると \(\dfrac{dx}{dt}=\alpha x\)
初期条件:1920年の人口 \(x=5600\) (万人) 比例定数:\(\alpha=0.0124\) とします。
\(t\) を経過年として、初めの10年間の人口を数値計算します。
\(x(0)=56\) 、 \(x(1)=56+0.0124\times 56=56.7\)
(赤の計算結果を確かめましょう)
時刻(年) | 1920 | 1921 | 1922 | 1923 | 1924 | 1925 | 1926 | 1927 | 1928 | 1929 | 1930 | \(\cdots\) |
人口(百万) | 56 | 56.7 | 57.4 | 58.1 | 58.8 | 59.6 | 60.3 | 61.0 | 61.8 | 62.6 | 63.3 | \(\cdots\) |
②計算で予測した人口をプロットします。
実人口の調査は5年おきにしか行われていないため、実人口のグラフは階段状になっています。
実人口の推移
時刻(年) | 1920 | 1925 | 1930 | 1935 | 1940 | 1945 | 1950 | 1955 | 1960 | 1965 | 1970 |
人口(百万) | 56 | 60 | 64 | 69 | 73 | 72 | 83 | 89 | 93 | 98 | 104 |
1940年代に実人口の停滞が見られるのは、第二次世界大戦が原因と考えられます。
理論式は \(x(t)=56e^{0.0124t}\) であり、計算予測とほぼ一致します。
Ⅴ 万有引力による惑星の軌道
太陽と惑星の間には、万有引力が働いています。
地球の重さを \(m_e\) 、太陽の重さを \(m_s\) 、距離を \(r\) 、万有引力定数を \(G=6.67259\times 10^{-11}\:[\rm N^2m^2/kg^2]\) として
\(F=\dfrac{Gm_em_s}{r^2}\)
地球の位置を \((x,y)\) 、太陽の位置を \((0,0)\) とすると、 \(\overrightarrow{F}=-\dfrac{Gm_em_s}{r^2}(x,y)\left. \right/ \sqrt{x^2+y^2}\)
引力は成分ごとに計算すればよいので、 \(F_x=-\dfrac{Gm_em_s}{\left(x^2+y^2\right)^{\frac{3}{2}}}x\:\:,\:\:
F_y=-\dfrac{Gm_em_s}{\left(x^2+y^2\right)^{\frac{3}{2}}}y\)
\(F_x=m_ea_x\:\:,\:\:F_y=m_ea_y\) なので、 \(a_x=-\dfrac{Gm_s}{\left(x^2+y^2\right)^{\frac{3}{2}}}x\:\:,\:\:
a_y=-\dfrac{Gm_s}{\left(x^2+y^2\right)^{\frac{3}{2}}}y\)
地球の軌道は、2階微分方程式
\(\dfrac{d^2x}{dt^2}=-\dfrac{Gm_s}{\left(x^2+y^2\right)^{\frac{3}{2}}}x\:\:,\:\:
\dfrac{d^2y}{dt^2}=-\dfrac{Gm_s}{\left(x^2+y^2\right)^{\frac{3}{2}}}y\) で表されます。
数値計算法を用いて、地球の軌道を算出しましょう。
太陽と地球の距離を、149600000000メートルとします。
天文学では、\(1.496\times 10^{11}\) メートル、エクセルでの表記は、\(1.496\rm E+11\) となります。
初期条件として、
地球の質量:\(m_e=5.974\rm E+24\:[\rm kg]\) 太陽の質量:\(m_e=1.989\rm E+30\:[\rm kg]\)
地球の初期座標:\((x,y)=(1.496\rm E+11\:,\:\:\:0\: )\)
地球の初期速度:\((v_x\:,v_y)=(0\:,\:\:2.978\rm E+04)\)
この条件のもとで、42日間の地球軌道を計算すると次の表のようになります。
(天文学的数値なので、電卓では計算しないでください)
時刻 [s] |
時刻 [h] |
時刻 [日] |
地球のx座標 | 地球のy座標 | ||||
位置 | 速度 | 加速度 | 位置 | 速度 | 加速度 | |||
x | vx | ax | y | vy | ay | |||
0 | 0 | 1.496E+11 | 0.000E+00 | -5.930E-03 | 0.000E+00 | 2.980E+04 | 0.000E+00 | |
600000 | 167 | 7 | 1.475E+11 | -3.558E+03 | -5.971E-03 | 1.788E+10 | 2.980E+04 | -7.240E-04 |
1200000 | 333 | 14 | 1.432E+11 | -7.141E+03 | -5.920E-03 | 3.550E+10 | 2.937E+04 | -1.468E-03 |
1800000 | 500 | 21 | 1.368E+11 | -1.069E+04 | -5.770E-03 | 5.259E+10 | 2.849E+04 | -2.219E-03 |
2400000 | 667 | 28 | 1.283E+11 | -1.415E+04 | -5.516E-03 | 6.888E+10 | 2.715E+04 | -2.962E-03 |
3000000 | 833 | 35 | 1.178E+11 | -1.746E+04 | -5.156E-03 | 8.411E+10 | 2.538E+04 | -3.681E-03 |
3600000 | 1000 | 42 | 1.055E+11 | -2.056E+04 | -4.690E-03 | 9.801E+10 | 2.317E+04 | -4.359E-03 |
原理を理解するため、簡略化した数値で試してみましょう。
地球の質量 \(m_e=1\:\rm[kg]\)
太陽の質量 \(m_s=10000\:[\rm kg]\)
地球の初期座標 \(x=10\:[\rm m]\) \(y=0\:[\rm m]\)
地球の初期速度 \(v_x=0\:[\rm m/s]\) \(v_y=1\:[\rm m/s]\)
万有引力定数を \(G=0.001\:[\rm N\cdot m^2/kg^2]\) とすると
\(a_x=-\dfrac{10}{\left(x^2+y^2\right)^{\frac{3}{2}}}x
\:\:,\:\:
a_y=-\dfrac{10}{\left(x^2+y^2\right)^{\frac{3}{2}}}y\)
\(t=0\) のとき \((x^2+y^2)^{\frac{3}{2}}=1000\)
\(a_x=-\dfrac{10}{1000}10=-0.1\) \(a_y=-\dfrac{10}{1000}\cdot0=0\)
\(t=1\) のとき
\(v_x=0+(-0.1)\times 1=-0.1\) \(\therefore x=10+(-0.1)\times 1=9.9\)
\(v_y=1+0\times 1=1\) \(\therefore y=0+1\times 1=1\)
\((x^2+y^2)^{\frac{3}{2}}=(9.9^2+1^2)^\frac{3}{2}=985.087\) より
\(a_x=-\dfrac{10}{985.087}\cdot 9.9=-0.1\) \(a_y=-\dfrac{10}{985.087}\cdot1=-0.01\)
\(t=2\) のとき
\(v_x=-0.1+(-0.1)\times 1=-0.2\) \(\therefore x=9.9+(-0.2)\times 1=9.7\)
\(v_y=1+(-0.01)\times 1=0.99\) \(\therefore y=1+0.99\times 1=1.99\)
\((x^2+y^2)^{\frac{3}{2}}=(9.7^2+1.99^2)^\frac{3}{2}=970.745\) より
\(a_x=-\dfrac{10}{970.745}\cdot 9.7=-0.10\) \(a_y=-\dfrac{10}{970.745}\cdot1.99=-0.02\)
次のような表が出来あがります。
時刻 | 地球のx座標 | 地球のy座標 | ||||
秒 | 位置 | 速度 | 加速度 | 位置 | 速度 | 加速度 |
[s] | x | vx | ax | y | vy | ay |
0 | 10 | 0 | -0.1 | 0 | 1 | 0 |
1 | 9.90 | -0.10 | -0.10 | 1.00 | 1.00 | -0.01 |
2 | 9.70 | -0.20 | -0.10 | 1.99 | 0.99 | -0.02 |
電卓で計算するのさえ大変ですが、もしあなたがアポロ13に乗船し全電源が停止していたなら、このような計算を手計算で行う必要があるでしょう。
エクセルのワークシートで続きを計算すると、地球の軌道は次の図のようになります。
このように地球の質量が1kgで太陽との距離が10mという世界があったなら、地球は太陽の周り約1分間で1周することになります。
[続く]