オイラー法を用いて空気抵抗を考慮した雨粒の落下に関する運動方程式を解いてみる

　さて、数学や物理学を学んだことのある方に、微分方程式を解いたことのある方も少なくないだろう。しかし世の中には我々が微分積分学で学んだ
解法では解くことが困難であったり、はたまた手計算では解けない微分方程式も存在する。そのようなものは実際存在するのか疑問に思う方もいるだろう。
では、以下の微分方程式をみてみよう。



図1. 解けない常微分方程式の例



この微分方程式は一見解けそうに見えるが実は初等函数では解けない（特殊函数になる）。一度手計算をしてみてはいかがだろうか。
以上のような微分方程式は一般に非線型の一階常微分方程式と呼ばれるもので、この微分方程式を解く場合、まず与えられた常微分方程式の斉次方程式を解くことから始まる。
しかしこの微分方程式は、その斉次方程式の段階から解けないのである。（解けるかもしれないが）

しかしこのような微分方程式をオイラー法や修正オイラー法と言った数値計算法で解くこともできる。ただしこれらの手法はあくまで『数値的に』解く手法であり、
手法によっては誤差が生じてしまうため精度が低いものもある。（オイラー法の場合テイラー展開の1次の項までしか取らないため精度が悪い。）


今回は、先ほど上に示した微分方程式ではないが、空気抵抗を考慮した雨粒の落下に関する運動方程式を、オイラー法を利用して解いてみようと思う。
ここでまず初めに運動方程式を示しておこう。



図2. 空気抵抗を考慮した雨粒の落下運動に関する運動方程式



ちなみに、この微分方程式は手計算で普通に解けます。今回の目的は、手計算で求めた場合と、オイラー法を用いた場合における誤差の検証を目標とする。
（運動方程式のグラフは『雨粒　運動方程式　グラフ』とググると出ると思います。申し訳ありません。）

次に数値計算に使用したプログラムを示そう。まずは1つ目のプログラムを以下に示す。



図3. プログラムその1



雨粒の重さは4[mg](4e-6[kg])、雨粒の初速度v_dは0[m/s]（時刻t=0[s]の時の速度）、重力加速度g=9.8[kg/s]、円周率π=3.1415、空気抵抗η=1.8e-5[Pa･s(N･s/㎡)]、
雨粒の半径r=1[mm](1e-3[m]、直径2mm）とする。また、b=6πηrとする。（物理学に関するサイトを参考にさせて頂きました。）

次にこのプログラムの出力結果を示そう。と言いたいところだが、データが多いためこれを記録したデータをこちらからダウンロードしていただきたい。（DLパスワードはcompile1です。）
プログラムを見ていただけるとお分かりだろうが、0.01[s]刻みで数値計算を行っており、またt=1.00[s]で計算を終了しているため、
これだけではまだグラフの概形を断定できない。（出力結果をエクセルにコピーしグラフを作成していただけるとなおよい。）そこで次に、少々プログラムを改良しt=50[s]の時まで
数値計算を行うプログラムを再び作成した。そのプログラムも以下に示そう。



図4. プログラムその2



同様にこのプログラムの出力結果も示したいところだが、前述の通りデータが多いため記録したデータをこちらからダウンロードしていただきたい。（DLパスワードはcompile2です。）
そして、この出力結果をgnuplotを用いて作成したグラフを以下に示すので見てみよう。



図5. プログラム2の出力結果から得られた運動方程式のグラフの概形



なお、縦軸はv、横軸はtである。このグラフを見ると最初の数[s]付近では急激に増加しているが時間が経過するにつれて速度の増加率が減少し、やがて一定の速度（終端速度）に収束する
ことが見て取れる。しかしこのグラフをよく見ると誤差がちゃんと見られます。



図6. 運動方程式の解



この運動方程式の解はこのようになる。ここでm,g,bの値は一定で、右辺における変数はtのみである。時刻tが経過するにつれて、exp(-bt/m)の値は常に減少していく（単調減少）ことは明らかである。
また、v(t)をtで微分すれば加速度v'(t)（a(t)でもいいが今回はv'(t）とする）が求まるので以下にその式を示す。



図7. 加速度v'(t)



gの値は一定で、また前述の通りexp(-bt/m)は単調減少であるので、加速度v'(t)は単調減少である。この事から時刻tが経過するにつれて微分方程式の解v(t)、すなわち雨粒の速度は増加はするが
加速度は徐々に落ちていき終いには一定の速度（終端速度）に収束していくことがわかるだろう。
ここで着目する点は、雨粒が一定の速度に収束すること、ではなく、『時刻tが経過するにつれてexp(-bt/m)は常に減少していくこと』である。
1-exp(-bt/m)は、t=0のとき1-1=0となり、十分に時刻が過ぎれば1-0=1に収束する。ここでexp(-bt/m)は単調減少でかつexp(-bt/m)≦1なので、1-exp(-bt/m)が単調増加であることは明らかであるが、
図5を見ると時刻tが45〜50[s]付近で速度が一度増加し、その後減少し再び一定の値に収束していることが見て取れる。（プログラム2の出力結果も参照していただけるとなお良い。）
これよりv(t)には極大値が存在することになるが、これは加速度v'(t)が0でない（すなわちv(t)が単調増加である）ことに矛盾する。（exp(-bt/m)はtがどのような値でも0にならないことは明らか。）
また同様に2[s]付近を見るとv(t)の変曲点が見られるように思えるが、これはv(t)の2階導関数v''(t)（躍度、または加加速度という。j(t)でもいい）が0でない（すなわち変曲点が無い）ことに矛盾する。

以上のことからオイラー法で運動方程式を解いた場合、理論値と比較するとやや誤差が生じるという結果に終わった。
プログラミング言語を用いて数値計算を行う際、int型やdouble型ではオーバーフローを起こす恐れがあるため、long int型やlong double型等を用いて再度数値計算を行う必要があると考えられる。
また、オイラー法では精度が低く大きな誤差が出る恐れがあるため修正オイラー法やルンゲ・クッタ法を用いた数値計算で再度検証を行う必要があると思われる。


〜あとがき〜

さて、プログラムを書いたり書き直したりgnuplotの使い方に悪戦苦闘したりグラフを出力したりで時間を多く食われ、とりあえず楽しんだもののこちらに結果（レポートとは違うか？）を実装するのに
更に時間を食われました。（編集が終わった頃は4時半を回っていた…）自分は数学が好き（元々勉学をする上で一番ましだった科目が数学だった）で、最近は高校数学のように答えだけを導く
だけでなく、初心に帰って理学部数学科で学ぶような数学を根本から学び直したり（大学院は数学専攻へ進学希望なので）していますが、現在通っている大学で専攻している学問で学びそうな数学、
即ちコンピュータを用いた数値計算（一応、情報工学を扱っているので）を楽しんでみることも、ひょっとしたら面白いのではないかと思いやってみたらやはり面白い…　数学は面白く、素晴らしい学問
だと思います。視覚的に見えない場所でも人間の生活をカバーしているのですから。無論、純粋に問題を解くことや、定理の証明とにらみ合いすることもありますが、このように趣向を変えてみるのも
たまには良いものだなと。実際、2年後期に数値計算は扱ったのだけれども、少々慣れない微分方程式で…面白かったけれども悪戦苦闘した記憶がある。というより今回は物理学の分野と絡めていた
けれども、2年後期に学んだものは生物学と絡めたものだったから、やや理解に苦しむものがあったのかもしれない。生物学は面白いけれどあまり勉強したことない（苦手生物いて抵抗ある）が、
微分方程式は数学、物理学だけでなく生物学や医学等の様々な分野で応用されているから、今後はもっと研究のし甲斐のある分野に手を付けてみたいと思う。そもそも自分が数学を専攻する
大学院へ進んだところで、何を専攻するかすら明確でないが。今回は近似計算をしたまでだけれども、あと1ヶ月で3年次に進級し、この学年で研究室への配属も決定する。所属希望の研究室は
ほぼ明確であるが、何を研究するかが明確でないため、今後、更に微分方程式に関する書籍を読む必要があるかもしれない。まず、数学専攻に進む場合は基礎からのやり直しが前提なのだがね。
さて、あとがきも長くなってしまうと読みづらくなりそうなので、今回はここまで。Studyページに関する次回の更新は未定です（蹴

※4月12日に一部修正しました。


2015/03/02 by プル23

戻る