ベルヌーイの微分方程式｜線形化の手順・特異解・リッカチとの関係

ベルヌーイの微分方程式は $y^{\prime }+P(x)y=Q(x)y^n$ という形をしています。右辺の $y^n$ のせいで非線形ですが、$z=y^{1-n}$ と置くだけで 1 階線形に変わります。

非線形の方程式が、従属変数を取り替えるだけで線形に移る数少ない例です。線形化という発想の原型にあたります。

以下では形の確認から始め、置換が $y^{1-n}$ になる理由を 3 つの角度から見たうえで、11 個の例を 1 つずつ解いていきます。

ベルヌーイ方程式の形

次の形の微分方程式をベルヌーイの微分方程式と呼びます。

$$\frac{dy}{dx}+P(x)y=Q(x)y^n$$

$P$ と $Q$ は $x$ の関数、$n$ は定数です。左辺は $y$ について 1 次ですが、右辺に $y^n$ があるため全体としては線形になりません。

$n$ は整数でなくてかまいません。$1/2$ や $2/3$ や $-1$ でもよく、応用では分数のべきがよく現れます。

n = 0 と n = 1 を除く理由

$n$ が $0$ と $1$ の場合は、置換を持ち出すまでもなく解けます。

置換 $z=y^{1-n}$ を見ると事情がはっきりします。$n=1$ では $z=y^0=1$ となって情報が消えるので、そもそも使えません。

$n=0$ では $z=y$ で、置換が何もしていないことになります。以後は $n\ne 0,1$ とします。

直観 1: 相対成長率で見る

両辺を $y$ で割ってみます。$y\ne 0$ のところでは次のように書き直せます。

$$\frac{y^{\prime }}{y}+P(x)=Q(x)y^{n-1}$$

左辺の $y^{\prime }/y$ は相対成長率で、$(\log y{)}^{\prime }$ に等しい量です。この式は、相対成長率が $y^{n-1}$ に応じて変わると読めます。

$n=1$ なら右辺から $y$ が消え、相対成長率が $y$ によらなくなります。指数関数的な増減になるのはこのためです。$n\ne 1$ のときだけ、大きさ自身が成長率を左右する構造が生まれます。

ロジスティック方程式が $n=2$ になるのも同じ理由です。個体数が増えるほど相対成長率が下がるという飽和の効果が、$y^{n-1}=y$ の項に入っています。

直観 2: 残るべきが 1-n になる

置換が $z=y^{1-n}$ になる理由は、割り算の結果を見れば分かります。両辺を $y^n$ で割ります。

$$y^{-n}\frac{dy}{dx}+P(x)y^{1-n}=Q(x)$$

$P$ の項に残った $y$ のべきがちょうど $1-n$ です。ここを 1 つの文字で置けば線形になる、というだけの話になります。

しかも第 1 項の $y^{-n}{y}^{\prime }$ は $y^{1-n}$ の微分の定数倍です。$(y^{1-n}{)}^{\prime }=(1-n)y^{-n}{y}^{\prime }$ なので、同じ $z$ で 2 つの項が同時に片づきます。

つまり $1-n$ という指数は選んだのではなく、割り算のあとに現れたものを拾っただけです。覚えるべきは指数ではなく手順のほうになります。

直観 3: 座標を取り替える

3 つ目の見方は、従属変数の座標変換です。$y$ という目盛りで見ると曲がっている方程式が、$z=y^{1-n}$ という目盛りでは真っすぐになります。

写像 $y\mapsto y^{1-n}$ は $y>0$ の範囲で全単射なので、その範囲では解の情報が失われません。$y=0$ と $y<0$ の扱いだけ別に考える必要があります。

うまい座標を選ぶと線形になる、という発想は他の非線形方程式にも広く使われます。最後のほうの節でいくつか挙げます。

線形化の手順

ここまでの内容を計算として書き下ろします。$z=y^{1-n}$ と置くと、微分は次のようになります。

$$\frac{dz}{dx}=(1-n)y^{-n}\frac{dy}{dx}$$

両辺を $y^n$ で割った式に代入します。

$$\frac{1}{1-n}\frac{dz}{dx}+P(x)z=Q(x)$$

$1-n$ を掛けて整理すると、$z$ についての 1 階線形微分方程式が得られます。

$$\frac{dz}{dx}+(1-n)P(x)z=(1-n)Q(x)$$

積分因子で解く

線形になったあとは機械的です。

積分因子を掛けた式は次のようになります。

$$\frac{d}{dx}\left[\mu (x)z\right]=(1-n)\mu (x)Q(x)$$

両辺を積分して $z$ を求め、最後に $y$ へ戻します。手順全体をまとめると 4 段です。

置換の対応表

よく出る $n$ について、置換がどうなるかを並べておきます。

$n>1$ なら $1-n<0$ なので、置換は $y$ の負べき、つまり逆数の形になります。$n<1$ なら正べきです。

もっともよく現れるのは $n=2$ の $z=1/y$ です。ロジスティック方程式もリッカチ方程式もここに落ちます。

例: 最も簡単な場合

簡単な例から始めます。

$$\frac{dy}{dx}+y=y^2$$

$n=2$ なので $z=y^{-1}$ と置きます。両辺を $y^2$ で割ると $y^{-2}{y}^{\prime }+y^{-1}=1$ です。

$z^{\prime }=-y^{-2}{y}^{\prime }$ なので、代入すると $-z^{\prime }+z=1$、整理して $z^{\prime }-z=-1$ になります。積分因子は $\mu =e^{-x}$ で、$(e^{-x}z{)}^{\prime }=-e^{-x}$ となります。

積分すると $e^{-x}z=e^{-x}+C$、つまり $z=1+C{e}^x$ です。$y=1/z$ に戻して次を得ます。

$$y=\frac{1}{1+C{e}^x}$$

失われる解

前節の答えには抜けがあります。$y\equiv 0$ も元の方程式の解ですが、どの $C$ を選んでも $1/(1+C{e}^x)$ はこれになりません。

原因は割り算です。両辺を $y^2$ で割った瞬間に、$y=0$ の場合が議論から外れています。

$n>0$ であれば $y\equiv 0$ はつねに $y^{\prime }+Py=Q{y}^n$ の解です。左辺も右辺も $0$ になるからです。

したがって一般解を書くときは、置換から出てくる解の族に加えて $y\equiv 0$ を特異解として添える必要があります。割り算をしたらその値を確かめる、という原則がここでも効いています。

例: 変数係数の場合

係数が $x$ に依存しても手順は変わりません。

$$\frac{dy}{dx}+\frac{y}{x}=x{y}^2$$

$n=2$ なので $z=y^{-1}$ です。両辺を $y^2$ で割ると $y^{-2}{y}^{\prime }+x^{-1}{y}^{-1}=x$ となり、代入して $-z^{\prime }+z/x=x$、つまり $z^{\prime }-z/x=-x$ になります。

積分因子は $\mu =\exp \left(-\int dx/x\right)=1/x$ です。両辺に掛けると $(z/x{)}^{\prime }=-1$ となり、積分して $z/x=-x+C$、すなわち $z=Cx-x^2$ です。

$y=1/z$ に戻して次を得ます。

$$y=\frac{1}{x(C-x)}$$

$n=2>0$ なので、これに特異解 $y\equiv 0$ が加わります。

例: n = 3 の場合

$n$ が 3 のときを見ます。力学系でもよく現れる形です。

$$\frac{dy}{dx}=y-y^3$$

$y^{\prime }-y=-y^3$ と書き直すと $n=3$ で、$z=y^{-2}$ と置きます。両辺を $y^3$ で割ると $y^{-3}{y}^{\prime }-y^{-2}=-1$ です。

$z^{\prime }=-2y^{-3}{y}^{\prime }$ なので $-\frac{1}{2}{z}^{\prime }-z=-1$、整理して $z^{\prime }+2z=2$ になります。積分因子 $e^{2x}$ を掛けて積分すると $z=1+C{e}^{-2x}$ です。

$y^{-2}=1+C{e}^{-2x}$ から次が得られます。

$$y=\pm \frac{1}{\sqrt{1+C{e}^{-2x}}}$$

三つの平衡解

前節の解の振る舞いを見ます。$x\to \infty$ で $e^{-2x}\to 0$ なので $y\to \pm 1$ となります。

<svg viewBox="0 0 700 306" xmlns="http://www.w3.org/2000/svg" style="max-width:100%;height:auto;display:block;margin:0 auto"><line x1="90" y1="160" x2="620" y2="160" stroke="#8a8a8a" stroke-width="1" stroke-dasharray="5 4"/><line x1="90" y1="74" x2="620" y2="74" stroke="#8a8a8a" stroke-width="1" stroke-dasharray="5 4"/><line x1="90" y1="246" x2="620" y2="246" stroke="#8a8a8a" stroke-width="1" stroke-dasharray="5 4"/><path d="M 110 132 C 250 132, 300 84, 600 76" fill="none" stroke="#1b81e0" stroke-width="1.5"/><path d="M 110 34 C 260 34, 320 68, 600 72" fill="none" stroke="#1b81e0" stroke-width="1.5"/><path d="M 110 188 C 250 188, 300 236, 600 244" fill="none" stroke="#1b81e0" stroke-width="1.5"/><path d="M 110 288 C 260 288, 320 252, 600 248" fill="none" stroke="#1b81e0" stroke-width="1.5"/><text x="634" y="78" font-size="12" fill="#1f1f1f">y = 1</text><text x="634" y="164" font-size="12" fill="#1f1f1f">y = 0</text><text x="634" y="250" font-size="12" fill="#1f1f1f">y = -1</text><text x="300" y="152" font-size="12" fill="#8a8a8a">不安定な平衡解（解はここから離れる）</text><text x="120" y="66" font-size="12" fill="#8a8a8a">安定</text><text x="120" y="266" font-size="12" fill="#8a8a8a">安定</text></svg>

$y\equiv 0$ も解ですが、まわりの解がすべて離れていく不安定な平衡解です。図では $y>0$ から出た解が $+1$ へ、$y<0$ から出た解が $-1$ へ向かいます。

割り算で失われた特異解が、力学系としては重要な役割を持っています。特異解を書き落とすと、解の全体像を取り違えることになります。

例: n = 1/2 の場合

指数が分数でも同じ手順です。

$$\frac{dy}{dx}+y=\sqrt{y}$$

$n=1/2$ なので $z=y^{1/2}$ です。両辺を $y^{1/2}$ で割ると $y^{-1/2}{y}^{\prime }+y^{1/2}=1$ になります。

$z^{\prime }=\frac{1}{2}{y}^{-1/2}{y}^{\prime }$ なので $2z^{\prime }+z=1$、整理して $z^{\prime }+\frac{1}{2}z=\frac{1}{2}$ です。積分因子 $e^{x/2}$ を掛けて積分すると $z=1+C{e}^{-x/2}$ になります。

$y=z^2$ に戻して次を得ます。

$$y={\left(1+C{e}^{-x/2}\right)}^2$$

$n=1/2>0$ なので、ここでも $y\equiv 0$ が特異解として加わります。

例: n = -1 の場合

$n$ が負の場合を見ます。事情が少し変わります。

$$\frac{dy}{dx}+y=\frac{1}{y}$$

$n=-1$ なので $z=y^{1-(-1)}=y^2$ です。両辺に $y$ を掛けると $y{y}^{\prime }+y^2=1$ になります。

$z^{\prime }=2y{y}^{\prime }$ なので $\frac{1}{2}{z}^{\prime }+z=1$、整理して $z^{\prime }+2z=2$ です。積分因子 $e^{2x}$ を掛けて積分すると $z=1+C{e}^{-2x}$ になります。

$y=\sqrt{z}$ に戻して次を得ます。

$$y=\sqrt{1+C{e}^{-2x}}$$

ここでは $y\equiv 0$ が解になりません。元の方程式が $y=0$ で定義されていないからです。

n の符号で変わること

直前の 2 例が示すように、$n$ の符号で状況が変わります。

$0<n<1$ のときはさらに注意が要ります。$y^n$ が $y=0$ で微分できないので、解の一意性が崩れることがあります。

この点はあとで $y^{\prime }=\sqrt{y}$ の例で確かめます。

ロジスティック方程式

応用でもっともよく現れるのがロジスティック方程式です。

$$\frac{dy}{dx}=ry\left(1-\frac{y}{K}\right)$$

$r$ は増殖率、$K$ は環境収容力を表します。展開すると $y^{\prime }-ry=-\frac{r}{K}{y}^2$ となり、$n=2$ のベルヌーイ方程式だと分かります。

$z=y^{-1}$ と置きます。両辺を $y^2$ で割ると $y^{-2}{y}^{\prime }-r{y}^{-1}=-r/K$ で、代入すると $-z^{\prime }-rz=-r/K$ です。整理して次になります。

$$\frac{dz}{dx}+rz=\frac{r}{K}$$

ロジスティック曲線の導出

前節の線形方程式を解きます。積分因子は $\mu =e^{rx}$ です。

$(e^{rx}z{)}^{\prime }=\frac{r}{K}{e}^{rx}$ を積分して $e^{rx}z=\frac{1}{K}{e}^{rx}+C$、すなわち $z=\frac{1}{K}+C{e}^{-rx}$ になります。

$y=1/z$ に戻すと、よく知られた形が出ます。

$$y=\frac{K}{1+CK{e}^{-rx}}$$

初期値 $y(0)=y_0$ を入れると $CK=(K-y_0)/y_0$ です。$x\to \infty$ で $y\to K$ となり、S 字型の曲線を描きます。

<svg viewBox="0 0 700 300" xmlns="http://www.w3.org/2000/svg" style="max-width:100%;height:auto;display:block;margin:0 auto"><line x1="90" y1="252" x2="620" y2="252" stroke="#1f1f1f" stroke-width="1"/><line x1="90" y1="30" x2="90" y2="262" stroke="#1f1f1f" stroke-width="1"/><line x1="90" y1="68" x2="620" y2="68" stroke="#8a8a8a" stroke-width="1" stroke-dasharray="5 4"/><path d="M 108 249 C 232 249, 254 168, 332 160 C 410 152, 452 76, 600 71" fill="none" stroke="#1b81e0" stroke-width="2"/><path d="M 108 36 C 250 36, 300 62, 600 66" fill="none" stroke="#1b81e0" stroke-width="1.5"/><circle cx="332" cy="160" r="4" fill="#1f1f1f"/><text x="634" y="72" font-size="12" fill="#1f1f1f">y = K</text><text x="634" y="256" font-size="12" fill="#1f1f1f">y = 0</text><text x="348" y="156" font-size="12" fill="#8a8a8a">変曲点（y = K/2）</text><text x="150" y="60" font-size="12" fill="#8a8a8a">K より大きい初期値からは減少して K へ</text><text x="150" y="238" font-size="12" fill="#8a8a8a">小さい初期値からは S 字を描いて K へ</text></svg>

$y\equiv 0$ と $y\equiv K$ の 2 つが平衡解です。前者は不安定、後者は安定で、$0$ と $K$ のあいだから出た解はすべて $K$ に近づきます。

変曲点は $y=K/2$ のところにあります。増加の勢いがもっとも強い時点で、そこから先は飽和に向かいます。

von Bertalanffy の成長式

分数のべきが現れる実例として、魚や動物の成長を表す von Bertalanffy の式を挙げます。

$$\frac{dw}{dt}=\alpha w^{2/3}-\beta w$$

体重 $w$ の増加が、体表面積に比例する同化の項 $w^{2/3}$ と、体重に比例する異化の項 $w$ の差で決まる、という考え方です。

$w^{\prime }+\beta w=\alpha w^{2/3}$ と書けば $n=2/3$ で、$z=w^{1/3}$ と置きます。両辺を $w^{2/3}$ で割ると $w^{-2/3}{w}^{\prime }+\beta w^{1/3}=\alpha$ です。

$z^{\prime }=\frac{1}{3}{w}^{-2/3}{w}^{\prime }$ なので $3z^{\prime }+\beta z=\alpha$、整理して $z^{\prime }+\frac{\beta }{3}z=\frac{\alpha }{3}$ になります。

成長曲線の形

前節を解きます。積分因子 $e^{\beta t/3}$ を掛けて積分すると $z=\frac{\alpha }{\beta }+C{e}^{-\beta t/3}$ です。

$w=z^3$ に戻して次を得ます。

$$w={\left(\frac{\alpha }{\beta }+C{e}^{-\beta t/3}\right)}^3$$

$t=0$ で $w=0$ とすると $C=-\alpha /\beta$ になり、次の形に整います。

$$w={\left(\frac{\alpha }{\beta }\right)}^3{\left(1-e^{-\beta t/3}\right)}^3$$

$t\to \infty$ で $w\to (\alpha /\beta {)}^3$ に収束します。これが漸近的な最大体重にあたり、水産資源の解析で使われる標準的な成長曲線です。

リッカチ方程式との関係

リッカチ方程式は次の形をしています。

$$\frac{dy}{dx}=q_0(x)+q_1(x)y+q_2(x)y^2$$

$q_0=0$ なら $y^{\prime }-q_{1}y=q_2{y}^2$ となり、$n=2$ のベルヌーイ方程式そのものです。定数項の有無が両者を分けています。

特殊解 $y_1$ が分かっているとき、$u=y-y_1$ と置きます。差を計算すると次のようになります。

$$\frac{du}{dx}-\left(q_1+2q_2{y}_1\right)u=q_2{u}^2$$

$n=2$ のベルヌーイ方程式です。特殊解 1 つで、解けない方程式が解ける方程式に変わります。

リッカチの例

実際にやってみます。

$$\frac{dy}{dx}=y^2-\frac{2}{x^2}$$

$y_1=a/x$ の形を試すと $-a/x^2=(a^2-2)/x^2$ となり、$a^2+a-2=0$ から $a=1$ または $a=-2$ が出ます。$y_1=1/x$ を使います。

$u=y-1/x$ と置くと $u^{\prime }=y^{\prime }+1/x^2=y^2-1/x^2$ で、$y=u+1/x$ を代入すると $u^{\prime }=u^2+2u/x$ になります。

$u^{\prime }-\frac{2}{x}u=u^2$ はベルヌーイ方程式です。$z=u^{-1}$ と置くと $z^{\prime }+\frac{2}{x}z=-1$ となり、積分因子 $x^2$ を掛けて積分すると $x^{2}z=-x^3/3+C$ です。

$u=1/z$ から $u=3x^2/(A-x^3)$（$A=3C$）となり、次を得ます。

$$y=\frac{1}{x}+\frac{3x^2}{A-x^3}$$

答えの検算

前節の答えを確かめておきます。$A=0$ とすると $y=1/x-3/x=-2/x$ です。

これは最初に見つけたもう 1 つの特殊解に一致します。$y^{\prime }=2/x^2$ と $y^2-2/x^2=4/x^2-2/x^2=2/x^2$ が等しいので、たしかに解になっています。

$A\to \infty$ の極限では第 2 項が消えて $y=1/x$、つまり出発点にした特殊解に戻ります。

2 つの特殊解が、一般解の両端として現れました。特殊解を 1 つ知っているだけで解の族全体が復元される、という構造が見てとれます。

一意性が壊れる例

$0<n<1$ では解の一意性が保証されません。もっとも簡単な例で確かめます。

$$\frac{dy}{dx}=\sqrt{y}$$

$P=0$、$Q=1$、$n=1/2$ のベルヌーイ方程式です。変数分離でも解けて $2\sqrt{y}=x+C$、すなわち $y=(x+C{)}^2/4$ が出ます。

さらに $y\equiv 0$ も解です。ここまでは今までと同じですが、この 2 つを貼り合わせた関数もまた解になります。

<svg viewBox="0 0 700 272" xmlns="http://www.w3.org/2000/svg" style="max-width:100%;height:auto;display:block;margin:0 auto"><line x1="90" y1="216" x2="640" y2="216" stroke="#1f1f1f" stroke-width="1"/><line x1="120" y1="40" x2="120" y2="230" stroke="#1f1f1f" stroke-width="1"/><line x1="120" y1="216" x2="620" y2="216" stroke="#1b81e0" stroke-width="2.5"/><path d="M 240 216 Q 384 216 512 66" fill="none" stroke="#1b81e0" stroke-width="1.5"/><path d="M 340 216 Q 466 216 580 76" fill="none" stroke="#1b81e0" stroke-width="1.5"/><path d="M 440 216 Q 546 216 636 100" fill="none" stroke="#1b81e0" stroke-width="1.5"/><circle cx="120" cy="216" r="4.5" fill="#1f1f1f"/><circle cx="240" cy="216" r="3.5" fill="#1f1f1f"/><circle cx="340" cy="216" r="3.5" fill="#1f1f1f"/><circle cx="440" cy="216" r="3.5" fill="#1f1f1f"/><text x="112" y="238" font-size="12" fill="#1f1f1f" text-anchor="middle">O</text><text x="240" y="238" font-size="12" fill="#8a8a8a" text-anchor="middle">a1</text><text x="340" y="238" font-size="12" fill="#8a8a8a" text-anchor="middle">a2</text><text x="440" y="238" font-size="12" fill="#8a8a8a" text-anchor="middle">a3</text><text x="150" y="196" font-size="12" fill="#1b81e0">y = 0 のまま進み、好きな点で立ち上がれる</text><text x="360" y="56" font-size="12" fill="#8a8a8a">どれも同じ初期条件 y(0) = 0 を満たす</text></svg>

$x\le a$ で $y=0$、$x>a$ で $y=(x-a{)}^2/4$ と定めると、$x=a$ で微分可能につながり全体で方程式を満たします。$a$ はいくらでも選べるので、原点を通る解は無限にあります。

原因は $\sqrt{y}$ が $y=0$ でリプシッツ連続でないことです。$n\ge 1$ や $n\le 0$ ではこの問題は起きません。

一般解の公式が書けることと、初期値問題の解が 1 つに定まることは別だ、という点を示す例になっています。

ベルヌーイでない例

形が似ていても置換が効かない方程式があります。対比のために挙げておきます。

判定は簡単です。$y$ について 1 次の項と $n$ 乗の項の 2 種類だけで書けているかを見ます。

3 種類以上の $y$ のべきが現れたら、ベルヌーイの手法は使えません。

線形化という考え方

ベルヌーイ方程式の教訓は、うまい従属変数を選べば非線形が線形になる場合がある、という点です。同じ発想は他でも使われます。

どれも解の空間を線形にするのではなく、方程式を線形にする座標を探す、という発想です。ベルヌーイ方程式はその最初の成功例にあたります。

もっとも、こうした変換が見つかる方程式はごく一部です。だからこそ、見つかったときの威力が際立ちます。

歴史

この方程式は 17 世紀末に提出され、すぐに解かれました。

微分積分学が生まれて間もない時期の出来事です。新しい道具でどこまで解けるかを探っていた時代の、代表的な成果の 1 つになります。

Bernoulli 兄弟は互いに問題を出し合う関係にありました。同じ時期に懸垂線や最速降下線の問題も扱われています。

参考文献