１階線形は積微分と部分積分で

このワークシートはMath by Codeの一部です。 今回も、方程式のタイプに応じた解法となる、基本的な視点をふやそう。微分方程式で式が単純なのは、もちろん微分が1階（1回）だけで、式も線形のものだ。だから、１階線形のタイプについて、さらに視点をふやしてみよう。 yのｎ階微分をy(n)と書き、Akがxの関数、ΣAk*y(k)=0のとき、線形同次。 Akが定数なら、特に、定数線形同次というね。線形だとしても、さらに、＝０ではなく、=B(x)のとき、それぞれ非同次という。

微分方程式y'=-y+e^2xは１階線形

１．1階線形

＜定数変化法＞ 定数変化法は同次方程式の一般解の定数を変数として変化させて、非同次方程式の解を求める方法。たとえば、dy/dx-y=xはdy/dx=x+yとなり、y/xの式ではないので、同次ではないね。そこでx=0としておきあとで復活させる。・dy/dx=y。1/y dy=1 dxで変数分離すると、∫1/y dy=∫1 dxによって、log|y|=x＋c 　 Cを定数とする一般解はy=Ce^x・このCをｘの関数として、c(x)とおき、もとの式に代入しよう。y=c(x)e^xとおく。このyを最初の式に入れると、dy/dx-y=(c(x))'e^x+c(x)(e^x)'-c(x)(e^x)=(c(x))'e^x=x。最後の等式を変形して、dc(x)/dx=x/e^x=x・e^-x。c(x)=∫e^-xx dxとなる。 部分積分する。∫Gf = FG- ∫fGとして、 f,G=e^-x,xとするとF,g=-e^-x,1 から、Fg＝-e^-xとなり積分はe^-xだから、c(x)=x・(-e^-x)-e^-x+D。つまり非同次の一般解はy=(x・(-e^-x)-e^-x+D)e^xy=-x-1+De^x ＜1階線形微分方程式＞ 上と同様にしてdy/dx+p(x)y=q(x)についても、いったんq(x)=0として、∫1/ydy=-∫p(x)dx=r(x)とおくと、(r(x))'=-p(x)となる。 log|y|=r(x)+C、Cを定数とする一般解がy=Ce^r(x)。C=c(x)として、y=c(x)e^r(x)積の微分(FG)'= fG+Fgから、 dy/dx+p(x)y=(c'(x)e^r(x)+c(x)(e^r(x))')+p(x)c(x)e^r(x)=q(x) r(x)'=-p(x)によって相殺される部分がある。 c'(x)e^r(x)=q(x) だから、c(x)=∫ q(x)・e^-r(x)dx+D。つまり、 y=e^r(x)(∫ q(x)・e^-r(x)dx+D)　　　r(x)=-∫p(x)dxとするとき（別解） dy/dx+p(x)y=q(x) r=∫p(x)dxとおき、両辺にe^rをかけると、 dy/dx e^r + p(x) y e^r = q(x) e^r 左辺は積の微分法則を逆算すると、y' (e^r) + y (p(x) e^r )=(ye^r)' だから、(ye^r)' =q(x) e^r つまり、ye^r = ∫q(x) e^r dx+c r=∫p(x)dxとするときy=1/e^r(∫ q(x)・e^rdx+c)　　　（例） y'+y=e^2xp(x)=1,q(x)=e^2xから、∫p(x)dx=∫1dx=xだから、 y=1/e^x(∫e^2xe^x+c)=1/e^x(∫e^3x+c)=1/e^x(1/3e^3x+c) y=1/3e^2x+c/e^x

微分方程式y'+y = xy^2はベルヌーイ型

微分方程式y'+y = xy^3はベルヌーイ型

２．ベルヌーイ型

ベルヌーイ型は線形1階に似ているけれど、xだけの項の代わりにｘだけの項にyⁿをかけた形 dy/dx+p(x)y=q(x)yⁿ の微分方程式だ。ｚ=y^1-n=y/yⁿと変数変換すると、1/(1-n)dz/dx+p(x)z=q(x)と、１階線形に直せる。（理由）合成関数の微分で、dz/dx=dz/dy dy/dx= (y^1-n) dy/dx =(1-n)/yⁿ dy/dx dy/dx= 1/(1-n) dz/dx yⁿ と変形できる。もとの式に代入すると、 1/(1-n) dz/dx yⁿ+p(x)y=q(x)yⁿ 両辺をyⁿで割る。 1/(1-n) dz/dx +p(x)ｚ=q(x) これで、ｚの一階線形に変換できたね。（例） y'+y = xy² 1-n=-1から、z=y^-1とおくと - z' + z =x から、 z' - z =-x から、 p(x)=-1, q(x)=-x r=∫p(x)dx=∫(-1)dx=-xとおくと、１階線形と部分積分で(integral(e^-x)=-e^-x,(-x)'=-1）から、 z=1/e^r(∫ q(x)・e^rdx+c)=1/e^-x(∫ (-x)・e^-xdx+c)=e^x(∫ (-x)・e^-xdx+c)=e^x{(-x)(-e^-x)-∫(-1)(-e^-x)+D)} ｚ=x+1+e^xD=y^-1 y=1/(x+1+e^xD) （例） y'+y = xy³ 1-n=-2から、z=y^-2とおくと (-1/2) z' + z =x から、 z' -2 z =-2x から、 p(x)=-2, q(x)=-2x r=∫p(x)dx=∫(-2)dx=-2xとおくと、１階線形, 部分積分で(integral(e^-2x)=-1/2e^-2x,(-2x)'=-2）から、 z=1/e^r(∫ q(x)・e^rdx+c)=1/e^-2x(∫ (-2x)・e^-2xdx+c)=e^2x{(-2x)(-1/2e^-2x)-∫(-2)(-1/2e^-2x)+D)} =e^2x{xe^-2x+1/2e^-2x+D)}= x + 1/2 + De^2x=y^-2 y²(x+1/2+e^2xD)=1

ロジスティックモデル

３．ロジスティック

指数関数のグラフは時間とともに急激に増大してしまう人口爆発予測の曲線として使われたりする。でも、実際は人口が無限大ということはなく、どこかで頭打ちがあるだろう。それも表せるのが ロジスティックモデルだ。 y'- Ay=- By²(A,Bは正の定数） これは、ベルヌーイ型でもあるね。 y'-Ay = -By² 1-n=-1から、z=y^-1とおくと - z' -A z =-B から、 z' + Az =B から、 p(x)=A, q(x)=B r=∫Adx=Axとおくと、１階線形と部分積分で(integral(e^Ax)=1/Ae^Ax,(B)'=0）から、 z=1/e^r(∫ q(x)・e^rdx+c)=1/e^Ax(∫ (B)・e^-Axdx+c)=1/e^Ax(∫ (B)・e^-Axdx+c) =1/e^Ax{(B)(1/Ae^Ax)-∫(0)(1/Ae^-Ax)+D)}=1/e^Ax{(B/A)e^Ax+D} z=(B/A) + D/e^Ax=y^-1 y=1/(B/A+D/e^Ax)となるね。 これは、ｘが時間、ｙが人口で B=0なら指数型人口増加のモデルy=(1/D)e^Ax(Dは正）なる。微分方程式はy'=Ayで、人口増加のスピードdy/dxは、そのときの人口y(x)に比例するというマルサスの法則にあてはまるね。 Bが正なら、制御項（- By²）が効いていて単調に変化して、極限値A/Bに収束する。 y(0)がA/Bより少ないときはDが正で単調増加し、y(0)がA/Bより多いときはDが負で単調減少する。

４．実装

質問：１階線形微分方程式をsympyで解くにはどうしたらよいでしょうか。 yが関数であることをy=Function('y')で宣言します。 y'+y=e^2xは、y'をDerivative(y(x),x)、yをy(x)で、e^2xをexp(2*x)とします。等式なので、eq=Eq(Derivative(y(x),x) +y(x), exp(2*x)と等式自体を変数にしましょう。これをy(x)の陽関数としての微分方程式を解くために dsolve(eq,y(x))としましょう。 y=1/3e^2x+c/e^x が求められるはずです。 #[IN]========================== from sympy.abc import * from sympy import * y = Function('y') eq = Eq(Derivative(y(x), x)+y(x), exp(2*x)) dsolve(eq, y(x)) #[OUT]=========================

だから、もちろん、ベルヌーイ型でも同じようにできるでしょう。なれてきたら、1行でかけますね。ｙをy(x)として書くことだけ気をつけよう。 y'+y = xy² y=1/(x+1+e^xD) #[IN]======================= dsolve( Eq(Derivative(y(x), x)+y(x), x*y(x)**2), y(x)) #[OUT]=========================