備忘録：非同次常微分方程式の解き方～未定係数法～

忘れてしまったのでぐぐってみたんだけど、いまいちピンとこれるとこがなかったので
自分にとってより具体的に理解できるようにブログにログを残してみるウェブ日記。


たとえば速度の1乗に比例した空気抵抗を受ける落下運動の式

mv'=-rv+mg

こいつは1階なので特に意識しなくても
mdv/dt=-rv+mgを適切に移項すれば

(-m/r)dv/dt=v-mg/r
dv/(v-mg/r)=(-r/m)dt
ln(v-mg/r)=-rt/m+C
v-mg/r=Cexp(-rt/m)
v=mg/r+Cexp(-rt/m)

初期条件t=0でv=0
v=0=mg/r+C
C=-mg/r
v=mg/r(1-exp(-rt/m))

が得られるんですよ。


でも、たとえば水か空気の抵抗を受けながら振動するバネとかだったら2階になるので
mx''=-kx-rx'+mg
こいつを解くときはさっきのようにはいかないんすよ。

この式は非同次というもので
mx''+rx'+kx=mg
の右辺がゼロじゃないんす。

まずは右辺をゼロにして、同次の微分方程式
mx''+rx'+kx=0

を解いてから非同次を解くっちゅうことになります。
この例における同次方程式は、重力によるつりあいをあらかじめ自然長に織り込み済みの式っちゅう具体的な状況に相当してます。

オーソドックスな微分方程式の解は、関数の何階微分が自分自身に比例とかそういう構造だけに
何階やっても変えれない指数関数exp(λt)を当てにすることが多いです。
先ほどの空気抵抗を受けた落下の式も、結果的に指数関数になったのはもはや必然っすね。


なので、とりあえずx=exp(λt)を微分方程式にぶち込んでみて
(mλ^2+rλ+k)exp(λt)=0
を満たすλを求めちまえっちゅう、このλのn次式を特性方程式と呼びます。


n階の微分方程式の特性方程式の解λは一般にn個ありますが
時々重解が現れるので注意です。


さて、λが2次方程式の解の公式でλ=(-r±√(r^2-4mk))/(2m)と求まったとしましょう。
λ1=(-r+√(r^2-4mk))/(2m)、λ2=(-r-√(r^2-4mk))/(2m)
とします。

exp(λ1t)もexp(λ2t)も関数xの条件を満たしているので
それぞれに適当な比を持たせても関数xは解だろうっちゅう推測の下、積分定数C1とC2を用いて
x=C1exp(λ1t)+C2exp(λ2t)
これが同次のときの一般解になります。



では非同次、つまり重力による垂れパイをxの自然長としてまだ考慮に入れてない式
mx''+rx'+kx=mg
だったらどうなのか

これを解くには未定係数法と定数変化法の2つがあります。
方法は違いますが同じ結果を導けます。


今回は未定係数法をやりましょう。

ただし、
まずは2階の微分方程式に入る前に、
腕試しのつもりで、空気抵抗を受けた落下の運動(1階の微分方程式)が、
さっきと同じ結果になることを確かめましょう。

mdv/dt=-rv+mgを移項して
mv'+rv=mg
のmgを外して同次にします。
mv'+rv=0

v=exp(λt)とおいて、式にぶち込みます。
(mλ+r)exp(λt)=0
なので、λ=-r/mです。

v=Cexp(-rt/m)　これが同次の一般解になりますが
次はいよいよ非同次の式

mv'+rv=mg
を解きましょう。


未定係数法では、
非同次の一般解yは
同次の解y0と非同次オンリーの解y1の線形結合になるので
適当な係数をかけて
その係数を求めるっちゅうのが未定係数法という方法になります。


非同次の右辺がある特定の形であれば
非同次オンリーの解y1は以下の表のようになることが知られています。

とにかく、右辺の指数の肩にある係数と特性方程式の解が一致したらとりあえず変数t掛け乗せしたれってこってす
重解だったらさらにtかけたれっちゅう。

mv'+rv=mgの右辺は表によるところのn次多項式×exp(αt)の、0次多項式・α=0の場合に相当し、
特性方程式の解λとαが不一致なので
非同次の解v1は適当な係数Aを設定してv1=Aとおいて

v1=Aとv1'=0を微分方程式にぶち込んで
0・m+rA=mg
よって、A=mg/rとなり、
一般解v=v0+v1=Cexp(-rt/m)+mg/rと、初期条件を決める前の式までちゃんと一致したことが確認できたと思います。
初期条件を求めるところからは先ほどと同じです。


次回は2階非同次の微分方程式を実際に未定係数法で解いてみましょう。
水とか空気とか油でダンパされたバネの、減衰する振動っすね。



ミケーネーｺｰﾎｰ　トゥービーコンティニュー

にほんブログ村