定数変化法をやろう
1。定数変化法
このワークシートはMath by Codeの一部です。
前回まで、定数係数非同次の微分方程式を未定係数法で解きました。
今回は、定数変化法で解いてみましょう。
その前に、微分方程式の具体にいろいろ触れてきたので、
まとめもかねて、一般論からやってみよう。
・同次式 y"+p(x)y'+q(x)y=0
一般解は線形独立な基底y1,y2の一次結合c1y1+c2y2。
もし、y1=ky2のように基底が比例する場合は、従属という。
ロンスキー行列式W=で判定できる。
W≠0となるx1があれば独立。
・非同次式 y"+p(x)y'+q(x)y=r(x)
一般解はy=y0+ys
y"+p(x)y'+q(x)y=0の一般解をy0=c1y1+c2y2、 y"+p(x)y'+q(x)y=r(x)の特解をysとする。
特解ysはy0のc1,c2に特定の値を入れて y"+p(x)y'+q(x)y=r(x)となるもの。
<定数変化法で線形方程式を解く>
2階線形非同次微分方程式y"+p(x)y'+q(x)y=r(x) p(x),q(x),r(x)は連続関数、y1, y2を解の基底
特解ysはys(x)= -y1∫y2 r /W dx + y2∫y1 r / W dx (Wはロンスキー行列式)
(例)
y" -3y’+2y=2x+1
同次式の一般解がt2-3t+2=0の解がt=1,2からy1=e1x,y2=e2x
ロンスキー行列式W=e1x(e2x)'-(e1x)'e2x=2e3x-e3x=e3x は正。
特解は定数変化法でys= -y1∫y2 r /W dx + y2∫y1 r / W dx
=-e1x∫(e2x)(2x+1)/e3xdx+e2x∫(e1x)(2x+1)/e3xdx
=-e1x∫(2x+1)e-xdx+e2x∫(2x+1)e-2xdx
=-e1x{(2x+1)(-e-x)+∫2e-xdx }+e2x{(2x+1)(-1/2e-2x)-∫(-e-2x)dx (部分積分した)
=-e1x{(2x+1)(-e-x)+(-2e-x) }+e2x{(2x+1)(-1/2e-2x)-(1/2e-2x)}
=(2x+1)+2 +{(2x+1)(-1/2)-1/2} =2x+3-x-1=x+2
非同次の一般解はys+y=x+2+Ae1x+Be2x
(例)
y" +2y’+y=4ex
同次式の一般解が(t+1)2=0の解t=-1から、y1=e-x, y2=xe-x
ロンスキー行列式W=e-x(xe-x)'-(e-x)'xe-x=e-x(e-x-xe-x)-(-e-x)xe-x
=e-2x-xe-2x+xe-2x=e-2x
特解が定数変化法でys= -y1∫y2 r /W dx + y2∫y1 r / W dx
=-e-x∫(xe-x)(4ex)/e-2xdx+xe-x∫(e-x)(4ex)/e-2xdx
=-4e-x∫xe2xdx+4xe-x∫e2xdx
=-4e-x{(1/2xe2x)-1/2∫(e2x)dx} +4xe-x(1/2e2x) (部分積分した)
=-4e-x{(1/2xe2x)-1/2(1/2e2x)} +4xe-x(1/2e2x)
=-2xex +ex +2xex
=ex
非同次の一般解はys+y= ex+Ae-x+Bxe-x
2.定数でないときも定数変化法で
未定係数法でやるのと同じ一般解がでましたね。
でも部分積分の嵐で、計算が大変でした。
定数変化法は部分積分をはじめとして、公式を駆使することが必要になります。
定数係数ならば未定係数法でカンタンに解けるので、あえて使わなくてもよいでしょう。
力を発揮するのは、係数がxの関数のときです。
<オイラー・コーシー方程式>
オイラー・コーシー方程式x2y'' + axy'+ by=0
は微分の階数分だけxの次数を上げているのでyはべき関数と予想がつくね。
y=xmとおいて、y'=mx(m-1), y''=m(m-1) x(m-2)から、方程式に代入すると、
左辺=xm{m(m-1)+am+b}=0 両辺をxmで割ると、特性方程式m2+(a-1)m+b=0ができる。
・mが2実数解m1,m2を持つと、y=c1xm1+c2xm2
・mが重複解をもつのはm=(1-a)/2のとき、y1=xm、y2=lnx y1とすると線形独立になる。
y=c1y1+c2y2
・mが複素数解m+ni, m-niを持つとき、x^in=(e^lnx)^in=e^(in lnx)から
x^(m+ni)=x^m e^(in lnx)=x^m(cos(n lnx) + i sin (n lnx))
x^(m-ni)=x^m e^(-in lnx)=x^m(cos(n lnx) - i sin (n lnx))
これから、y1= x^m(cos(n lnx)), y2=x^m(sin(n lnx))
(例)
x2y'' -2.5xy' -2y=0
特性方程式m2+(a-1)m+b=m2-3.5m-2=0 の解はm=-0.5,4 一般解はy=c1x-0.5+c2x4
(例)
x2y'' -3xy' +4y=0
特性方程式m2+(a-1)m+b=m2-4m+4=0 の解はm=2 一般解はy=x2(c1+c2 ln x)
(例)
x2y'' +7xy' +13y=0
特性方程式t2+(a-1)t+b=t2+6t+13=0 の解はt=-3±2i 一般解はy=x-3(Acos(2 lnx) + Bsin(2 ln x))
<定数変化法で解く>
オイラー・コーシー方程式の非同次版にとりくもう。
x2y'' -2xy' +2y=x3lnx
m2+(a-1)m+b=m2-3m+2=0 の解はm=1,2 解の基底はy1=x, y2=x2 で、一般解はy=c1x1+c2x2
ここで、方程式の両辺をx2で割り標準形にすると、
y′′−2 y′/x+2 y/x2=x ln(x) となり、非同次項をf(x)=xln(x) とします。
ロンスキー行列式 W=x(x2)'-(x)'x2
=2x2-x2=x2
特解は定数変化法で
ys= -y1∫y2 f(x) /W dx + y2∫y1 f(x) / W dx
=-x∫x2xlnx/x2 dx+x2∫x xlnx/x2 dx
= -x∫x*lnx dx +x2∫1*lnx dx
= -x{(lnx *1/2 x2)- ∫(1/x*1/2 x2) dx} +x2{(lnx *x)- ∫(1/x*x) dx} (部分積分)
= -1/2x{(lnx *x2)- ∫ x dx} +x2{(lnx *x)- ∫1 dx}
= -1/2x(lnx * x2)- 1/2 x2 ) +x2(lnx *x- x )
= -1/2 lnx * x3+ 1/4 x3 +lnx *x3- x3
= 1/2 lnx * x3 -3/4 x3
一般解はy+ys= 1/2 lnx * x3 -3/4 x3 +c1x1+c2x2
2.実装
質問:オイラー・コーシー方程式を解析的に解くコードはどうかいたらよいですか。
<Sympy>
sympyでは、同次と同じように非同次が解けるはずです。
y = Function('y')(x)とかいておくと、yがxの関数であることを宣言しているので、
それ以降いちいち、yをy(x)と書く必要がないようです。
微分方程式はEq(左辺, 右辺 )とかきます。
これをy=の形の一般解を得るためにnsolve(微分方程式, y)とかけばよいね。
x2y'' -2.5xy' -2y=0
x2y'' -3xy' +4y=0
x2y'' +7xy' +13y=0
x2y'' -2xy' +2y=x3lnx
を解いてみよう。
#[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(b(x), c(x), f(x), x開始値, y開始値, y’開始値, x終了値, 間隔)
2階の常微分方程式 y′′+b(x)y′+c(x)y=f(x)が解けます。
(例)
SolveODE(x^2, 2x, 2x^2 + x, x(A), y(A), 0, 5, 0.1)
は,あらかじめ定義した A を出発点として2階の常微分方程式を解く.
常に結果を軌跡として返す.アルゴリズムは現在のところ ルンゲ=クッタ法 に基づく. |
