ちょぴん先生の数学部屋

数学の楽しさを、現役メーカーエンジニアが伝授するぞ!

微分方程式への招待 ~ニュートンの運動方程式~

今回の記事では、近代物理学の至宝、「ニュートン運動方程式」について紹介し、実際にそれを使っていくつかの物理現象を解明していきます。

 

それにあたっては「微分方程式を解く」という作業が必要になりますので、各物理現象を通して、その一端に触れていければと思います。

 

0. ニュートン運動方程式とは?

 

高校物理で習う話ではありますが、改めて紹介しておくと、

「物体の加速度は、その物体に加わっている力と比例する」というニュートンによって発見された物理法則になります。いわゆる「F=ma」ってやつですね。

 

加速度というのは、速度を時間で微分してできるものであり、その速度は物体の変位(位置)を時間で微分することで求まります。まとめると、加速度は物体の変位を時間で2回微分したものになります。

 

なので大学以降では、ニュートン運動方程式は以下のように書かれるのが普通です。

力とか加速度は「大きさ」と「向き」の両方を持った量なので、このようにベクトルの形で書かれます。

 

さて、この式は「Fの情報から、変位xを時間tの関数として求める」というタイプの方程式になっていて、求める関数xに微分が絡んでいます。このような方程式を微分方程式と呼びます。

 

ニュートン運動方程式に限らず、多くの物理法則は微分方程式の形で書かれることが一般的で、個々のシチュエーションに合わせてその微分方程式を解くことで、現象を詳しく分析するというのが、物理学の基本スタイルになります。

 

1. 具体例(Fがxによらないとき)

 

まずは、力Fが物体の位置xに依存しない一定値を取る場合について考えます。

1-1. 自由落下(空気抵抗が無視できる場合)

図のように、高さHの位置から質量mの物体を落とす状況を考えます。

座標軸を上向きにとったとき、物体にかかる力は下向きの重力-mgだけなので、運動方程式は下のようになります。

この重力は常に一定なので、tにもyにも依存していません。

なので、yを求めようと思ったら、素直に時間tで2回積分すればOKです。

すると、高校物理でお馴染みの式が求まります。

 

特に初速v0が0のとき、地上に落ちる時刻Tは以下の式で計算できることになります。

 

1-2.ボール投げ

 

図のように原点の位置から、質量mのボールをθの角度で初速v0で投げる状況を考えると、運動方程式は以下のようになります。

今回はx,yという2方向を考える必要があるので、運動方程式はこのようなベクトルの形になるわけです。

ボールにかかる力はy方向下向きの重力mgだけで、x方向には何も力がかかっていません。

 

ベクトルの形になっているとは言え、x方向についての考察とy方向の考察は独立しているので、まずは別個に考えてしまいましょう。

 

初速v0ですが、これをx方向とy方向に分解するとそれぞれ、v0cosθ, v0sinθとなることに注意して方程式を各成分ごとに解くと以下のようになります。

 

これでボールの座標がtの式で書けたので、tを消去してあげればボールの軌跡が求まります。

 

このようにボールの軌跡は「上凸の放物線」になります。といいますか、ボールを投げた時にできる図形が2次関数になるので、2次関数のグラフの事を「放物線」と呼ぶわけです。

 

さて、ボールを投げた結果、ボールがx=Lの位置で着地するとすると、Lは以下のように計算できます。

このようにLはsin2θに比例するので、飛距離Lを最大にしようと思ったら、2θ=90°⇔θ=45°の角度で投げればよいことになります。

 

2. 具体例(Fがxに依存するとき)

 

次に、力Fの式にxが含まれている場合を考えていきます。

 

2-1. ばね振動

図のように、ばね定数kのばねに質量mの重りをつけて、自然長からAだけ伸ばした状態で手を放す状況を考えると、運動方程式は以下のようになります。

 

座標は、自然長のときにx=0として、伸びる方向を正にとっています。このとき、重りに働く力は、伸びた方向と逆向きで、フックの法則からkxと書けます。

 

このように、力がxの関数になっている場合、先ほどまでのような単純な積分で解ける、というわけにはいきません。

 

さて、この運動方程式をよく見てみると、「xを2回微分しても(係数を除き)xの形のまま」という性質を持っていることが分かります。

 

微分をしても形が変わらない関数、といえば「指数関数」ですね。

 

ということで、xを以下のような指数関数の形で仮定してしまいます。

 

これを元の方程式に代入すると、以下のようにλが求まります。虚数になってしまいましたね・・・とりあえずそれは脇に置きます。

λが2つ求まった場合、xは一般的に以下のように、各λに対応した指数関数の定数倍の和の形(以後、「線形和」と呼ぶことにします)に書くことができます。

実際に元の方程式に代入すれば成り立ってることがすぐに分かるはずです。

 

初期条件は、t=0のときx=A, dx/dt=0なのでC1, C2が求まります。

これでxの式が求まったのですが、虚数が残ってて気持ち悪いですね。ここは、オイラーの公式数学界のKingとQueenは、愛で結ばれた・・~世界で一番美しい数式、オイラーの等式~ - ちょぴん先生の数学部屋 (hatenablog.com))を使って、三角関数の形に直してしまいましょう。

これでシンプルになりました。

xは振幅A、周期2π/ωの単振動の式になり、高校物理の公式通りの結果になります。

 

2-2. 抵抗があるばね振動

次に、2-1で考えたばねとおもりを、粘性νの液体に漬けて同じ動作を行う状況を考えます。すると、おもりにはkxだけでなく、速度dx/dtに比例した抵抗力が働くことになります。

以上を踏まえて運動方程式を立てると以下のようになります。

抵抗力の符号がマイナスになっている理由は、おもりの進行方向(つまり速度dx/dt)と逆向きに抵抗力が働くからです。

 

この方程式についても、xを指数関数の形で仮定すると考えやすくなります。2-1と同様にλを計算すると、以下のようになります。

ルートの中身が正か負かによって、λが実数か虚数かが変わってくるので、場合分けします。

 

[Ⅰ]ルートの中身が正の時

 

ν^2>4mkのときで、いわば「粘性が強すぎるとき」に相当します。このときλは2つの負の実数になるので、それぞれ-a,-bとおくと、xは指数関数の線形和でかけたので、以下のようになります。

初期条件からC1,C2を確定させてあげると、下のようになります。

このように、xは指数の肩が負の値をとる指数関数の形で書けて、tについて単調減少しt→∞でx=0に収束することが分かります。これは速度dx/dtについても同じです。

要するに、粘性が強すぎて、振動することなくおもりの動きはすぐに止まってしまうことになります。

 

[Ⅱ]ルートの中身が0の時

 

このときλは重解になるので、このままだと線形和が作れません(2回微分する微分方程式なので、積分定数が2つ必要なはずです)。ここは、天下り的ではあるのですが、

指数関数の「係数」をtの関数だとみなして、元の方程式に代入するとうまくいきます。

 

これでC1,C2と積分定数が2つにできました。

初期条件からこの2つを確定させると、以下のようになり、

[Ⅰ]の場合と同じくxは単調減少して0に収束する関数になります。このようなλが重解になるときの挙動を「臨界減衰」と呼び、想像通りに非常にレアな現象です。

 

[Ⅲ]ルートの中身が負の時

 

ν^2<4mkのときで、いわば「粘性が弱いとき」に相当します。このときλは共役な虚数になります。[Ⅰ]と同じように指数関数の線形和でかくと、xは以下のようになります。

指数関数の内、方が実数の物は共通しているので括れることになります。

 

初期条件からC1,C2を求めて、オイラーの公式虚数を消去すると、下のようになります

これは、振幅が指数関数的に小さくなりながら、振動する関数(減衰振動)になります。粘性が弱いがゆえに振動できる余力があるということです。

 

このように、粘性νの大きさによって、ばねは3通りの振る舞いをすることが分かります。

 

2-3. 単振り子

図のように長さLの紐の先に質量mのおもりをぶら下げた振り子を、振れ角θ0の状態から手を放す状況を考えます。

 

座標軸を、おもりの軌跡となる円周に沿って設定してあげると、運動方程式は以下のようになります。

おもりには重力と糸からの張力が働くわけですが、おもりは円周と垂直な方向には動かないので、今回はこの向きの運動方程式は無視しています。そうすると、円周に沿った向きの力は重力の成分mgsinθだけです。

 

今回の力はθの三角関数の形なので、2-1,2のばね振動の時のような「指数関数の形を仮定する」が使えません。

 

とはいえ、θが小さい値の時は近似が可能なので、近似解と厳密な解の両方を考えてみることにします。

 

[Ⅰ]近似解

 

θが小さい値の時は、sinθ≒θと近似できます。この近似によって、単振動と同じような式に帰着出来て、綺麗に解くことができます。

特に、この振り子の周期Tは下のようになり、

ひもの長さLだけに依存し、振れ幅θ0と質量mには依存しないことが分かります。これは、ガリレオが発見した「振り子の等時性の性質そのものです。

 

あくまで上記は近似解ではありますが、周期については、実際はθ=45°くらいの比較的大きい角度に対しても厳密解との誤差が4%程度に収まる、割と精度の良い近似になっています。

 

[Ⅱ]厳密解

 

前述したように「指数関数の形に仮定」の解法が使えないので、一工夫を加えます。

 

両辺に試しにdθ/dtをかけて見ましょう、すると、「合成関数の微分」の性質から、

と言った感じに、微分をひとまとめにできます。

このとき、この式は、「(tの関数)の微分=0」となっています。これは、微分する前の(tの関数)が、実はtに依存しない定数である、ことを意味しています。つまりこういうことです。

初期条件からCを決め、θの範囲を0°~90°に限定すると以下のようになります。

2乗を外すときに-をつけているのは、θが単調減少するためです。

 

ここからなのですが、残念ながらθをtの式で綺麗に解くことはできません。

 

なのでθを直接解くことはあきらめて、周期Tだけ求めることにします。

 

今、θが-Δθだけ変化するときの時間経過Δtは、微分方程式から

と分かります。なので、両辺をθ=0~θ0の範囲で全部足してしまえば、時間経過の合計は周期の1/4になるはずです。微小な量の足し算こそが積分なので、以下のようになります。

最後の式の積分は綺麗に解けないので、これ以上簡単にはなりません(実際に求めようと思ったら数値計算が必要)。

このように、厳密解では周期が振れ幅θ0にも依存していることが分かるので、振れ幅が大きすぎると「振り子の等時性」は成り立たなくなります。

 

この厳密解の周期について、

と置換し、さらに近似解での周期をT0とすると、厳密解と近似解の比は次のようになります。

この積分の形は「第1種完全楕円積分」と呼ばれるもので、数値計算しやすい形になっています。実際にpythonを使って数値計算すると下のようなグラフになり、

実際には振れ幅が45°の場合でも誤差が4%程度でかなり精度がよく、振れ幅90°の場合は誤差18%で流石に少し大きいですが、近似解は思ったよりは悪くない近似式だとわかります。

 

(参考)第1種完全楕円積分数値計算するpythonコード

 

2-4. 空気抵抗のある自由落下

最後に、冒頭の1-1に空気抵抗を追加した状況を考えます。初期状態はy=H, dy/dt=0とします。

 

空気抵抗の大きさは速度dy/dtに比例すると考えると、運動方程式は以下のようになります。

力に、yに依存するものとそうでないものが混在していますね。

 

まず、

方程式の解の1つになっていることがすぐに分かります。この解を使って、微分方程式の解が、

の形になっていると仮定し、このy0についての微分方程式を考えることにしましょう。

すると、

となって、y0に依存する項だけになったので、これで「指数関数の形を仮定」の解法に持ち込めるようになりました。y0を求めると下のようになるので、

本来求めたかったyは、

となります。

さらに初期条件からA1,A2を決めた最終的な答えは下のようになります。

そこから、ボールの速度は下のように求まります。

ボールの速さ|dy/dt|は単調増加して、t→∞でmg/γに収束することがわかります。つまり、ボールの速さはやがて一定値になってしまうという事です。

 

 

このように、ニュートン運動方程式を利用することで、様々な身近な物理現象を調べられる、というわけです。