定数係数2階線形非同次は未定係数法で
強制振動と共鳴１
１．定数係数2階線形非同次
このワークシートはMath by Codeの一部です。
前回まで、2階線形同次を判別式で分類して解きました。
今回は、同じ2階線形だけど、非同次の場合について調べてみよう。
＜非同次の特解＞
２階のa,bが定数の線形非同次微分方程式y"+ay'+by=f(x)の一般解をｙ、特殊解をysとする。
解法は特別解、特殊解(特解）が決め手になる。
なぜか？
微分の線形性から、非同次一般解yー非同次特殊解syは同次一般解y0になる。
y-ysは(y-ys)''+a(y-ys)'+b(y-ys)=y"+ay'+by-(ys''+ays'+bys)=f(x)-f(x)=0から、y"+ay'+by=0の解となる。
y- ys=y0となるから、y＝y0＋ys
つまり、非同次の一般解yは、同次一般解y0に特解ysを加えればよいのです。
＜未定係数法＞
f(x)が定数、べき関数、指数関数、三角関数なら、特殊解は予想しやすい。
そこで、未定係数をつけたy＝特殊解の候補を代入して、決定するのが未定係数法。
f(x)の種別による特解ysの求め方。
・定数f(x)=cなら、y''=y'=0を代入して0+a*0+by=cを解いて、ys=c/b
・1次式f(x)=px+qなら、y=Ax+Bを代入してa*(A)+b(Ax+B)=px+qの係数比較から
bA=p、aA+bB=qから、A=p/b, B=q/b^2-ap/b。ys=p/b x + (qb-ap)/b2
・指数関数f(x)=Becxなら、y=Aecxを代入して、 (Aecx)''+a(Aecx)'+bAecx=ecxから、ecxで割り、
(c2+ac +b)A=B A=B/(c2+ac +b) より、ys＝B/(c2+ac +b) ecx
・三角関数f(x)=Bcos(cx)またはBsin(cx)なら、y"+ay'+by=Beicxとおき、その特解を求める。
y=Aeicxの代入(-c2+iac +b)A=Bから、ys=B/(-c2+iac +b) eicxより、
f(x)=Bcos(cx)の特解はRe(ys),f(x)=Bsin(cx)の特解はIm(ys)
または、y=Asincx + Bcoscxを代入して係数比較する。
・f(x)がｎ次多項式なら、ysもｎ次多項式で未定係数を使います。
・他にf(x)が多項式、指数関数、三角関数の積ならば、ysも未定係数でそれらの積を代入します。
（例）
y" -3y’+2y=2x+1 特解がys=2/2 x +(1*2-(-3)*2)/(22)=x+2
同次式の一般解がt2-3t+2=0の解がt=1,2からy=Ae1x+Be2xだから、
非同次の一般解はys+y=x+2+Ae1x+Be2x
（例）
y" +2y’+y=4ex 特解がys=4/(12+2(1)+1) ex= ex
同次式の一般解が(t+1)2=0の解t=-1から、y=(Ax+B)e-1xだから、
非同次の一般解はys+y= ex+(Ax+B)e-x
（例）
y" -2y’+2y=5sin 2x
・yを複素関数とみなして、y" -2y’+2y=5ei2xとおく。
y=Aei2xを代入して、ei2xでわり(-4-4i+2)A=5 A=5/{(-2)(1+2i)}=-1/2(1-2i)から、
ys=Aei2x=-1/2(1-2i)(cos 2x + i sin2x) = -1/2{(cos 2x +2sin2x) + (-2cos2x+sin2x) i }
特解はIm(ys)=-1/2sin2x+cos2x
・または、y=Asin2x + Bcos2xとおき、代入する。
(Asin2x + Bcos2x)" -2(Asin2x + Bcos2x)’+2(Asin2x + Bcos2x)=5sin 2x
(-4Asin2x -4Bcos2x) -(4Acos2x - 4Bsin2x)+(2Asin2x + 2Bcos2x)=5sin 2x
sin2x(-4A+4B+2A-5) + cos2x(-4B -4A + 2B)=0
-2A+4B=5 -2B=4A B=-2A -10A=5 A=-1/2 , B=1 ys=-1/2sin2x + cos2x
特性方程式はt2-2t+2=0 (t-1)2=-1 t=1±i
同次式の一般解はy=ex(A cosx + B sinx)だから、
非同次の一般解はys+y=-1/2sin2x + cos2x+ex(A cosx + B sinx)
２．実装
質問：非同次2階定数係数の微分方程式を解析的に解くコードはどうかいたらよいですか。
＜Sympy＞
sympyでは、同次と同じように非同次が解けるはずです。
y = Function('y')(x)とかいておくと、ｙがｘの関数であることを宣言しているので、
それ以降いちいち、yをy(x)と書く必要がないようです。
微分方程式はEq(左辺, 右辺 )とかきます。
これをy=の形の一般解を得るためにnsolve(微分方程式, y)とかけばよいね。
「y" -3y’+2y=2x+1」
「y" +2y’+y=4ex 」
「y" -2y’+2y=5sin 2x」を解いてみよう。
#[IN]Python
#===========================
from sympy.abc import *
from sympy import *
y = Function('y')(x)
ddy=Derivative(y, x, 2)
dy=Derivative(y, x)
＜Geogebra＞
CASで非同次を解いてみよう。
SolveODE(微分方程式）とかきます。
３．実例
＜ばねの強制振動と共鳴＞
地震のショックをやわらげるダンパー、
橋をわたるときのゆれの増大、
いろんな場面で、ゆれを押さえたり増幅することがテーマになります。
これらは非同次微分方程式とつながりますよ。
前回、空気抵抗を入れた水平のばねを調べましたね。
空気抵抗は、ばねの先端のおもりの速さに比例して、時間とは無関係でした。
空気抵抗のかわりに時間とともに変化する力F(t)を作用させることを考えてみよう。
初期値は、ばねの自然の変位をx0、t=0の速度はv=0としましょう。
xを時間の関数としてF=mx''=-kx-F(t)。（ｍ、ｋは正）
つまり、mx''+kx=F(t)という非同次微分方程式ができるね。
でも、mがじゃまなので、
F(t)=mf cos(ωt) 。これは、振幅fで周波数ωで、質量mに比例する力だね。
k=m ω02。これは、単振動のときと同様で、ω02は正。
mx''+mω02x=mfcos(ωt)
x''+ω02x=fcos(ωt)
と簡単になる。
・一般解を求めよう。
x=Asinωt + Bcosωtとおき、代入する。
(Asinωt + Bcosωt)" +ω02(Asinωt + Bcosωt)=fcos(ωt)
(-ω2Asinωt -ω2Bcosωt) +(ω02Asinωt + ω02Bcosωt)=fcos(ωt)
sinωt(-ω2A+ω02A) + cosωt(-ω2B+ ω02B-f)=0
(ω02-ω2)A=0 (ω02-ω2)B=f
ω02-ω2≠0 なら、A=０、B=f/ (ω02-ω2)
特解はxs=f/(ω02-ω2)cosωt
特性方程式はt2+ω02=0 t=±i ω0
同次式の一般解はx=e0x(A cosω0t+ B sinω0t)=(A cosω0t + B sinω0t)だから、
非同次の一般解はxs+x=f/(ω02-ω2)cosωt + (A cosω0t + B sinω0t)
・初期条件を入れた特殊解は、
t=0で、v=0, x=x0。
x(0)=f/(ω02-ω2)cos0+ (A cos0 + B sin0)
=f/(ω02-ω2)+ A = x0 これから、A=x0- f/(ω02-ω2)
v(t)=dx/dt=(f/(ω02-ω2)cosωt + (A cosω0t + B sinω0t))'
=-fω/(ω02-ω2)sinωt +ω0(-A sinω0t + B cosω0t)
v(0) = -fω/(ω02-ω2)sin0 +ω0(-A sin0 + B cos0)
= ω0 B =0 だから、ω0 から, B=0
A,Bを代入して、
x(t)=f/(ω02-ω2)cosωt + (A cosω0t + B sinω0t)
＝f/(ω02-ω2)cosωt + (x0- f/(ω02-ω2)) cosω0t + 0 sinω0t)
＝ x0cosω0t + f(cosωt - cosω0t)/(ω02-ω2)
=
ここまでやってみて、何か気が付くことはないだろうか？
ω02-ω2≠0 であることは必要だけど、ω02とω2が限りなく近いときは意味があるよね。
質問：ω2がω02に限りなく近づくとどうなるでしょうか。
共鳴（レゾナンス）が起こります。
共鳴は、日常のコトバとしては、同じ考えに同調するという響きあうというよいイメージで使われたりしますね。
でも、物理的には、振動するものに外部から同じ振動数の力が加わったとき、波の山と山が重なり、増幅される現象です。
しかも、それが減衰せずに振幅が増大し続けると危険な結末になることもあるかもしれません。
だから、負のイメージですね。
みんなが吊り橋を同じリズムで渡ると橋がどんどん揺れるという怖い経験をしたことはありませんか。
数式でみると、特殊解のx(t)のfに乗じた分数部分は、ロピタルの定理により、ωの有理関数とみて微分してから、ωにω0を代入したものが極限値となる。
ωがω0に近づくとき、分数部分は、 が極限。
だから、ｘ(t)＝