微分方程式 解曲線の描画   ちょっと解説

このツールは2次元の微分方程式の解曲線を描画します。 割と簡単そうな方程式でも解析的に解くには結構骨が折れるし、 そもそも解析的に解くことができない場合もあります。 そこで、微分方程式をそのまま利用して近似解をwebで簡単に描画するツールを作りました。 近似は基本的に4次のrunge-kutta法で行っていますが、Euler法と修正Euler法も選択できます。 初期値は、第1象限にランダムに16点を選び、残りの3象限に反転コピーしています。 計64点を初期値として、一斉に描画が開始されます。

( 描画ページではLocal Storageを使用しています。脚注 [1] 参照。)

描画が開始されると現れる青色のラインは、初期値から通常の時間経過を経て描画されます。
しばらくすると、黄色いラインが現れます。こちらは初期値から時間を遡って描画されます。
青と黄色のラインの動きを観察することで、解曲線の性質を知ることができます。

それでは、描画ページのセレクトボックスに表示されるリストに従って、 ざっくりと説明します。

以下の説明は、森本光生氏の「パソコンによる微分方程式」からの引用を基にしています。

[01] 2次元1次線形自励系

\[ \frac{dx}{dt} = ax + by \] \[ \frac{dy}{dt} = cx + dy \]

係数入力ボックス $a, b, c, d$ に値を設定することで、 様々な解曲線が動的に描画されます。
試しに、以下の例の値を設定してみて下さい。 説明に沿った解曲線が描画されると思います。

ここで、$a, b, c, d$ は実数で 全ては0でないとして、 実係数の2次の正方行列\( A \) を導入します。

\[ A = \pmatrix{ a & b \cr c & d \cr } \]

$A$ の固有値を $\lambda_1, \lambda_2$ とします。
以下の説明で、「固有値に対応する直線」という表現が出てきますが、 これは各固有値の固有空間を意味しています。2行2列の正方行列の場合、 固有空間は一般に原点を通る直線になります。

  1. 固有値の1つが0の場合

    例)

    \[ \frac{dx}{dt} = x + y \] \[ \frac{dy}{dt} = 2x + 2y \]

    \[ A = \pmatrix{ 1 & 1 \cr 2 & 2 \cr } \]

    $\lambda_1 = 0$, $\lambda_2 = 3 > 0$より、 解曲線は $\lambda_2$ に対応する直線 $y = 2x$ に平行な半直線で、 その出発点は $\lambda_1$ に対応する直線$y = -x$ 上の各点であり、無限遠に向かう。

    $\lambda_1 = 0$, $\lambda_2 < 0$の場合、解曲線の向きは逆になる。

  2. 固有値が $\lambda_1 < 0 < \lambda_2$ の場合 (鞍状点)

    例)

    \[ \frac{dx}{dt} = -2x + 2y \] \[ \frac{dy}{dt} = -2x + 3y \]

    \[ A = \pmatrix{ -2 & 2 \cr -2 & 3 \cr } \]

    $\lambda_1 = -1 < 0 , \, \lambda_2 = 2 > 0$ より、 解曲線は $\lambda_2$ に対応する直線 $y = 2x$ に沿いながら無限遠から発し、原点方向に引き寄せられ、 あるところで$\lambda_1$ に対応する直線$2y = x$ に沿うように向きを変え、無限遠に向かう。

  3. $ 0 < \lambda_1 < \lambda_2$ の場合1 (不安定結節点 または 湧点)

    例)

    \[ \frac{dx}{dt} = x \] \[ \frac{dy}{dt} = 2y \]

    \[ A = \pmatrix{ 1 & 0 \cr 0 & 2 \cr } \]

    $\lambda_1 = 1 > 0 , \, \lambda_2 = 2 > 0$ より、 解曲線は $\lambda_1$ に対応する直線 $y = 0 \, (x軸)$ に沿って原点から湧き出し、
    $\lambda_2$ に対応する直線 $x = 0 \, (y軸)$ に沿って無限遠に向かう。

  4. $ 0 < \lambda_1 < \lambda_2$ の場合2 (不安定結節点 または 湧点)

    例)

    \[ \frac{dx}{dt} = x \] \[ \frac{dy}{dt} = 3x + 2y \]

    \[ A = \pmatrix{ 1 & 0 \cr 3 & 2 \cr } \]

    $\lambda_1 = 1 > 0 , \, \lambda_2 = 2 > 0$ より、 解曲線は $\lambda_1$ に対応する直線 $y = -3x$ に沿って原点から湧き出し、
    $\lambda_2$ に対応する直線 $x = 0 \, (y軸)$ に沿って無限遠に向かう。

  5. $ \lambda_1 < \lambda_2 < 0$ の場合 (安定結節点 または 沈点)

    例)

    \[ \frac{dx}{dt} = -4x - 3y \] \[ \frac{dy}{dt} = 2x -11y \]

    \[ A = \pmatrix{ -4 & -3 \cr 2 & -11 \cr } \]

    $\lambda_1 = -10 < 0 , \, \lambda_2 = -5 < 0$ より、 解曲線は $\lambda_1$ に対応する直線 $y = 2x$ に沿って無限遠から発し、
    $\lambda_2$ に対応する直線 $3y = x$ に沿って原点に収束する。

  6. 固有値 $ \lambda $ が重根で、行列Aが対角化可能な場合
    $ \lambda > 0$ のときは (不安定退化結節点 または 退化湧点)となる

    例)

    \[ \frac{dx}{dt} = 2x \] \[ \frac{dy}{dt} = 2y \]

    \[ A = \pmatrix{ 2 & 0 \cr 0 & 2 \cr } \]

    $\lambda_1 = \lambda_2 = 2 < 0 $ より、 解曲線は原点より発し、無限遠に至る半直線である。

    $\lambda < 0 $ であれば、解曲線は無限遠より発して原点に至る半直線になる(安定退化結節点 または 退化沈点)。

  7. 固有値 $ \lambda $ が重根で、行列Aが対角化不可能な場合
    $ \lambda > 0$ のときは (不安定退化結節点 または 退化湧点)となる

    例)

    \[ \frac{dx}{dt} = 2x + y \] \[ \frac{dy}{dt} = 2y \]

    \[ A = \pmatrix{ 2 & 1 \cr 0 & 2 \cr } \]

    固有値は重根で $\lambda = 2$ であるが、行列Aは対角化できない。 解曲線は原点より発し、固有値 $\lambda$ に対応する直線 $y = 0 \, (x軸)$ に沿って無限遠に至る。

    $\lambda < 0 $ であれば、解曲線は無限遠より発して固有値 $\lambda$ に対応する直線に沿って原点に至る。
    (安定退化結節点 または 退化沈点)。

  8. 固有値が純虚数の場合 (渦心点)

    例)

    \[ \frac{dx}{dt} = x + 2y \] \[ \frac{dy}{dt} = -x - y \]

    \[ A = \pmatrix{ 1 & 2 \cr -1 & -1 \cr } \]

    固有値は $\lambda = \pm i$ である。 解曲線は原点を中心とする楕円軌道(周期解)を描く。

  9. 固有値が実部と虚部を持つ場合 (渦状点)

    例)

    \[ \frac{dx}{dt} = -x - y \] \[ \frac{dy}{dt} = 2x - y \]

    \[ A = \pmatrix{ -1 & -1 \cr 2 & -1 \cr } \]

    固有値は $\lambda = -1 \pm \sqrt{2}i$ である。 固有値の実部が負の場合、解曲線は無限遠より発して原点の回りをらせん状に回りながら原点に至る、 安定渦状点である。

    固有値の実部が正の場合は、解曲線は原点より発して原点の回りをらせん状に回りながら無限遠に至る、 不安定渦状点である。

[02] 2次元2次線形自励系

\[ \frac{dx}{dt} = ax^2 + bxy + cy^2 \] \[ \frac{dy}{dt} = dx^2 + exy + fy^2 \]

係数入力ボックス $a, b, c, d, e, f$ に値を設定することで、 様々な解曲線が動的に描画されます。
ここで、$a, b, c, d, e, f$ は実数で 全ては0ではないとします。

森本氏の著書では、解析解を示され、解曲線の性質を詳しく解説されていますが、 ここでは、かなり端折った説明になっていることをご留意願います。

  1. 結節峠点型

    例)

    \[ \frac{dx}{dt} = x^2 \] \[ \frac{dy}{dt} = 2(x - y)^2 \]

    \[ \frac{dx}{dt} = 1 \cdot x^2 + 0 \cdot xy + 0 \cdot y^2 \] \[ \frac{dy}{dt} = 2 \cdot x^2 - 4 \cdot xy + 2 \cdot y^2 \] より、

    $ a = 1, \, b = 0, \, c = 0 $
    $ d = 2, \, e = -4, \, f = 2 $

    直線 $2y = x$ と $x = 0 \, (y軸)$ が2つの双曲型角領域を作っている。 また直線 $y = 2x$ に沿って4つの放物型角領域があり、そこでは解曲線がこの直線を挟んで、 原点に至る。第1象限にある放物型角領域は不安定であり、第3象限にある放物型角領域は安定である。

  2. 退化結節峠点型

    例)

    \[ \frac{dx}{dt} = 4x^2 \] \[ \frac{dy}{dt} = x^2 + 4y^2 \]

    \[ \frac{dx}{dt} = 4 \cdot x^2 + 0 \cdot xy + 0 \cdot y^2 \] \[ \frac{dy}{dt} = 2 \cdot x^2 + 0 \cdot xy + 4 \cdot y^2 \] より、

    $ a = 4, \, b = 0, \, c = 0 $
    $ d = 1, \, e = 0, \, f = 4 $

    直線 $2y = x$ と $x = 0 \, (y軸)$ が2つの双曲型角領域を作っている。 また直線 $2y = x$ に沿って2つの放物型角領域がある。 第1象限にある放物型角領域は不安定であり、第3象限にある放物型角領域は安定である。

  3. 双楕円型

    例)

    \[ \frac{dx}{dt} = -x^2 + y^2 \] \[ \frac{dy}{dt} = -2xy \]

    \[ \frac{dx}{dt} = -1 \cdot x^2 + 0 \cdot xy + 1 \cdot y^2 \] \[ \frac{dy}{dt} = 0 \cdot x^2 - 2 \cdot xy + 0 \cdot y^2 \] より、

    $ a = -1, \, b = 0, \, c = 1 $
    $ d = 0, \, e = -2, \, f = 0 $

    直線 $y = 0 \, (x軸)$ が平面を2つの楕円型角領域に分けている。 解軌道は原点で $x軸$ に接する円周で、原点より発し、原点に至る。

  4. 三峠型

    例)

    \[ \frac{dx}{dt} = x^2 - y^2 \] \[ \frac{dy}{dt} = -2xy \]

    \[ \frac{dx}{dt} = 1 \cdot x^2 + 0 \cdot xy - 1 \cdot y^2 \] \[ \frac{dy}{dt} = 0 \cdot x^2 - 2 \cdot xy + 0 \cdot y^2 \] より、

    $ a = 1, \, b = 0, \, c = -1 $
    $ d = 0, \, e = -2, \, f = 0 $

    双楕円型とは1つめの式の $x^2$ と $y^2$ の符号が変わっただけであるが、 解曲線はがらっと変わったものになる。 直線 $y = 0 \, (x軸)$ と$y = \pm\sqrt{3}$ によって 3つの尾根と谷が組み合わさった軌道になるため、3峠型と表現する。 3直線が平面を6つの双曲型の角領域に分離している。

  5. 結節点型 2

    例)

    \[ \frac{dx}{dt} = -x^2 \] \[ \frac{dy}{dt} = y^2 \]

    \[ \frac{dx}{dt} = 1 \cdot x^2 + 0 \cdot xy + 0 \cdot y^2 \] \[ \frac{dy}{dt} = 0 \cdot x^2 + 0 \cdot xy + 1 \cdot y^2 \] より、

    $ a = 1, \, b = 0, \, c = 0 $
    $ d = 0, \, e = 0, \, f = 1 $

    直線 $y = 0 \, (x軸) と x = 0 \, (y軸)$ が2つの双曲型角領域を作っている。 また直線 $y = -x$ に沿って4つの放物型角領域があり、そこでは解曲線がこの直線を挟んで、 原点に至る。第2象限にある放物型角領域は不安定であり、第4象限にある放物型角領域は安定である。

  6. 放物型

    例)

    \[ \frac{dx}{dt} = -x^2 \] \[ \frac{dy}{dt} = y^2 \]

    \[ \frac{dx}{dt} = -1 \cdot x^2 + 0 \cdot xy + 0 \cdot y^2 \] \[ \frac{dy}{dt} = 0 \cdot x^2 + 0 \cdot xy + 1 \cdot y^2 \] より、

    $ a = -1, \, b = 0, \, c = 0 $
    $ d = 0, \, e = 0, \, f = 1 $

    直線 $y = -x$ が平面を2つの放物型角領域に分けている。 直線 $y = 0 \, (x軸)$ が平面を2つの楕円型角領域に分けている。 解軌道は直線を挟んで、その向きが逆になる。

[03] 高次の同次自励系

高次の自励形以降は、係数の入力はできません。 事情はお察し頂けると思います。
ただし、 この方程式は係数が変えられたら面白いと思った場合は、最小限の係数を可変にしています。
最終的には、web上で自由に方程式を定義して、その解曲線を描画するページにしたいと思っていますが、 まだ先の話になりそうです。
今のところは森本氏の著書に紹介されている方程式を紹介したいと思います。

  1. 周期解を持つ同次自励形

    例)

    \[ \frac{dx}{dt} = -y^3 \] \[ \frac{dy}{dt} = x^3 \]

  2. 8つの放物型角領域を持つ同次自励形

    例)

    \[ \frac{dx}{dt} = x^3 \] \[ \frac{dy}{dt} = y^3 \]

    直線 $y = \pm x$ と $x, \, y 軸$ によって、平面は8つの放物型角領域に分けられる。

  3. 複合角領域を持つ同次自励形

    例)

    \[ \frac{dx}{dt} = x^2(y - x) \] \[ \frac{dy}{dt} = y^2(y - 2x) \]

    直線 $2y = (3 \pm \sqrt{5})x$ と $x, \, y 軸$ によって、 平面は8つの角領域に分けられる。双曲型と楕円型の角領域をそれぞれ2つ持ち、 放物型角領域を4つ持つ。

  4. Vinograd(ヴィノグラード)の自励形

    例)

    \[ \frac{dx}{dt} = x^2(y - x) + y^5 \] \[ \frac{dy}{dt} = y^2(y - 2x) \]

    前出の複合角領域を持つ同次自励形の $dx/dt$ の右辺に $y^5$ が加わったものだが、 $y = 0 \, (x軸)$ を境に点対称の形でムンクの叫びを連想させる楕円型角領域が現れる。

[04] 自励系の諸例

  1. x軸上に周期的に鞍状点と渦心点を持つ自励形

    例)

    \[ \frac{dx}{dt} = y \] \[ \frac{dy}{dt} = \sin x + y^2 \]

    この自励形は
    $(p\pi, 0) \quad p = \pm1, \pm2, \ldots $ に平衡点を持つ。
    $p$ が偶数のときは、鞍状点となり、奇数のときは、渦心点になります。 全体としては、金魚が並んだような解曲線が描かれます。奇数の渦心点が金魚の目になり、 偶数の鞍状点が尻尾になります。

    森本氏の丁寧な解説を思いっきり端折ります。 それぞれの平衡点の近傍で、$\sin x$を級数展開し、$1次$で近似します。
    $p$ が偶数のときは$p = 0$とすれば$x = 0$の近傍では$\sin x \approx x$と近似できます。 $y^2$ も$x$に比べて十分に小さいとすれば、 方程式は以下のように表現できます。 \[ \frac{dx}{dt} = y \] \[ \frac{dy}{dt} = x \] これを線形近似といいます。対応する2次の正方行列の固有値は$\pm 1$ になるので、平衡点は鞍状点になります。 [01] - 2 参照

    $p$ が奇数のときは、$p = \pi$とすれば$x = \pi$の近傍では$\sin x \approx -(x - \pi)$と近似でき、 $x^2$ も$x - \pi$に比べて十分に小さいとすれば、 方程式は以下のように表現できます。 \[ \frac{dx}{dt} = y \] \[ \frac{dy}{dt} = -(x - \pi) \] $x - \pi, -y$を新しい座標とすれば、 \[ \frac{dx}{dt} = -y \] \[ \frac{dy}{dt} = x \] 対応する2次の正方行列の固有値は$\pm i$ になるので、平衡点は渦心点になります。 [01] - 8 参照

    それぞれの平衡点の近傍では上記の近似が成立し、離れるに連れて、 元々の方程式の項の影響が現れ、全体として、金魚が並んだような解軌道になると考えられます。

  2. 円周上に8つの平衡点を持つ自励形

    例)

    \[ \frac{dx}{dt} = x^2y^2 - 1 \] \[ \frac{dy}{dt} = x^2 + y^2 - a \]

    この自励形は、円周 $x^2 + y^2 = a $ 上に8つの平衡点を持ちます。
    各象限に2つずつ計8個です。元々パラメーター$a$には固定値の 4 が与えられていましたが、 可変にしてみました。どのような解軌道が描かれるか、確認してみて下さい。 平衡点の位置が少し分かりづらいので、描画終了時に円を表示するようにしています。 デフォルトの描画範囲は 8 ですが、16 とか 32 にすると、 初期値をランダムに取っている関係で軌道は荒くなりますが、平衡点は確認し易くなると思います。

  3. 平面上に周期的な平衡点を持つ自励形

    例)

    \[ \frac{dx}{dt} = \cos x + \cos y \] \[ \frac{dy}{dt} = \sin x \sin y \]

    この自励形は $ \pmatrix{ 2p\pi \cr (2q + 1)\pi \cr } と \quad \pmatrix{ (2p + 1)\pi \cr 2q\pi \cr } に平衡点を持ちます(p, qは整数)。 $

    4-1と同様に、平衡点の近傍で1次近似を行うと、 1つ目のタイプの平衡点は双楕円型になります。
    2つ目のタイプの平衡点は3峠型になります。
    平面上で、これが繰り返されます。

  4. 極限軌道の例 1
    解軌道が極限円になる場合

    例)

    \[ \frac{dx}{dt} = ay + ax(1 - x^2 - y^2) \] \[ \frac{dy}{dt} = cx + dy(1 - x^2 - y^2) \]

    パラメーター$ a, b, c, d $ には以下の値が与えられています。

    $ a = 1, \, b = 1, \, c = -1 , \, d = 1 $

    極限軌道はPoincaré (ポアンカレ) によって導入されました。 この例では極限軌道は $x^2 + y^2 = 1$になり、極限円になります。 初期値が円の外側にある場合も内側にある場合も、円に螺旋状に巻きついて行きます。 この極限円は安定であるといえます。

    極限円については森本氏が結構詳しく解説されているのですが、私自身は理解できる状態に至っていません。 森本氏は "笠原皓司 微分方程式の基礎 朝倉書店 1982" を参考にすることを勧められています。

    さて、パラメーター$ a, b, c, d $ に以下の値をセットしてみて下さい。 解軌道の様相が一変します。

    $ a = -1, \, b = -1, \, c = 1 , \, d = -1 $

    今度も極限軌道は $x^2 + y^2 = 1$になりますが、 円の内側の解軌道は螺旋状に原点に向かい、 円の外側の解軌道は螺旋状に円から離れて無限大に向かいます。 従って、極限円は不安定です。

    極限軌道については、下記のサイトを参照してください。 次に紹介する予定の "von del Polの自励形" も含めて、詳しく解説されています。

    1. リミットサイクル

    "ポアンカレ 極限円 etc." で検索をかけていくと、"だまし絵で有名な Escher (エッシャー) が深く関係していることが見えてきます。

    2. ポワンカレの円板モデル

    こちらは幾何についての考察なのですが、何がどう絡んでいるのか...まさにカオスです。

  5. 極限軌道の例 2 von del Pol(ファン・デル・ポル)の自励形

    例)

    \[ \frac{dx}{dt} = y \] \[ \frac{dy}{dt} = -x + ay(1 - x^2) \]

    デフォルトではパラメーター$ a $ には 1 がセットされています。

    ファン・デル・ポルの方程式は様々なサイトで紹介されていますが、 [04] - 4 で参照先として紹介した リミットサイクル が分かりやすいと思います。 極限軌道は、円とはとてもいえないひしゃげた四角形に見えますが、 この方程式も、極限円の性質を備えているとのことです。 パラメーター $a$ を0.9, 0.8…と小さくして描画を実行すると、次第に円に近づいていくのが分かると思います。 このケースの解軌道も安定です。

    $ a = 0$ の場合、原点は2次元1次の線形自励形で紹介した渦心点 ( [01] - 8 ) になります。 更に、$a$ が負になるとどうなるでしょうか? 確認してみて下さい!

  6. 極限軌道の例 3 Lemniscate(レムニスケート)

    例)

    \[ \frac{dx}{dt} = a * y(1 + 2x^2 + 2y^2) \] \[ \frac{dy}{dt} = b * x(1 - 2x^2 - 2y^2) + c * y(1 + 2x^2 + 2y^2)((x^2 + y^2)^2 - (x^2 - y^2)) \]

    デフォルトではパラメーター$ a, b $ には 1 が、$c$ には -1 がセットされています。

    この自励形の平衡点は3つあり、原点 (鞍状点) と2つの安定な渦状点 ($ \pm 1 / \sqrt{2}, 0 $) である。

    レムニスケート と呼ばれる曲線 $ (x^2 + y^2)^2 - (x^2 - y^2) = 0 $ に、3つの極限軌道が存在する。 原点は平衡点で、$x < 0$ と$ x > 0$ の部分が方程式の解軌道になる。 森本氏によると、これらは極限円ではないとのこと。 これは2つの極限軌道が連結しているためと思われます。極限軌道自体は安定です。

    $ a = -1, b = -1, c = 1$ の場合はどうなるでしょうか? 確認してみて下さい!

    $ a = -1, b = 1, c = -1$ の場合はどうなるでしょうか? 3つの平衡点が原点 (渦状点) と2つの鞍状点 ($ \pm 1 / \sqrt{2}, 0 $) に変わります。 レムニスケートではなくなり、原点を中心に人の目のような軌道が現れます。
    $ a = 1, b = -1, c = 1$ では軌道の向きが反転します。

  7. 極限軌道の例 4 線形近似が純虚数の場合 1
    極限円が存在しない場合

    例)

    \[ \frac{dx}{dt} = y \] \[ \frac{dy}{dt} = -x - y^3 \]

    この自励形の平衡点は原点のみである。

    原点で線形化( 線形近似 )すると、式は以下のようになる。

    \[ \frac{dx}{dt} = y \] \[ \frac{dy}{dt} = -x \]

    固有値は $\lambda = \pm i$ になり、 解軌道は [1] - 8 で見たように渦心点になる。

    ところが、実際に元の非線形の方程式の解軌道をプロットしてみると、渦心点ではなく渦状点らしき軌道が現れる。

    この軌道についての森本氏の説明を要約すると次のようになる。

    原点Oを中心とする半径 $r$ の円周のy切片をA, Bとする。 この円周上の各点で $ ベクトル \pmatrix{ y \cr -x - y^3 \cr } は \quad 円周に接するベクトル \pmatrix{ y \cr -x \cr } と \quad 円の内側を向いているベクトル \pmatrix{ 0 \cr -y^3 \cr } の和になるので、 $
    円の内側を向いている。

    y軸上の点Aを通る解軌道は、Aより円周の内部に入り込み、 y軸上の点Bよりも原点Oに近いB'を横切る。 同じ議論をy軸上の点B'に置き換えると、 B'を通る解軌道は、B'より円周の内部に入り込み、 y軸上の点A'を横切る。

    こうして、Aを通る解軌道は原点の回りを螺旋状にぐるぐると回り閉じることなく 原点に収束する。

    ただし、原点近傍では軌道の動きが鈍くなり、 "刻み幅"を小さくして、"繰返し"を増やさないと 予測通りの軌道がプロットされなくなるので、確認の際は注意すること。

  8. 極限軌道の例 4 線形近似が純虚数の場合 2
    極限円が無限個現れる例

    例)

    \[ \frac{dx}{dt} = -y + 10x \, e^{(-1 / \sqrt{r})} \sin(1 / r)\] \[ \frac{dy}{dt} = x + 10y \, e^{(-1 / \sqrt{r})} \sin(1 / r)\] $\qquad$ where $\; r = \sqrt{x^2 + y^2} $

    $r$ は原点からの距離である。この方程式を(4.8.a)と呼ぶことにする。

    \[ \alpha(r) = \cases{ e^{(-1 / \sqrt{r})} \sin(1 / r) & ($r > 0$) \cr 0 & ($r = 0$) \cr } \]

    $\alpha(r)$ で定義される関数は原点に近づくに連れて減衰振動を繰り返す曲線である。

    添付の図が荒くて分かり難いが、 $r$のステップを小さくして、図を引き延ばせば、フラクタルであることが推察できる。

    $\alpha(r) = 0$ となる(グラフが横軸と交わる)位置に注目する。 後ほど触れるが、交点の座標は右側から原点に向かって$1 / \pi, 1 / 2\pi, 1 / 3\pi, 1 / 4\pi, 1 / 5\pi...$と続く。
    $1 / \pi ~ 1 / 2\pi$ の間では、$\alpha(r) < 0$ である。
    $1 / 2\pi ~ 1 / 3\pi$ の間では、$\alpha(r) > 0$ である。
    $1 / 3\pi ~ 1 / 4\pi$ の間では、$\alpha(r) < 0$ である。
    $1 / 4\pi ~ 1 / 5\pi$ の間では、$\alpha(r) > 0$ である。
    ...

    $\alpha(r) = 0$の場合、方程式(8.4.a)は右辺第2項がそれぞれ消えるので、線形自励形

    \[ \frac{dx}{dt} = -y \] \[ \frac{dy}{dt} = x \]

    と等しくなる。この方程式を(8.4.b)と呼ぶことにする。

    固有値は $\lambda = \pm i$ になり、 解軌道は [1] - 8 で見たように渦心点になる。

    従って、方程式(4.8.a)は原点を中心とする半径 $r = 1 / k\pi\quad(ただし、kは自然数)$ の円周上では、 線形化した方程式(4.8.b)と一致した閉じた解軌道(極限円)を持つ。

    $k$ は自然数であるから、閉じた解軌道は無限に存在することになるが、 軌道は半径 $1 / \pi$ の円周の内側に出現する。

    では、方程式(8.4.a)の初期値が円周 $1 / k\pi$ 上にない場合の解軌道は、どうなるだろうか? 方程式(8.4.a)から右辺第1項を消去した方程式(8.4.c)は

    \[ \frac{dx}{dt} = 10x \, e^{(-1 / \sqrt{r})} \sin(1 / r)\] \[ \frac{dy}{dt} = 10y \, e^{(-1 / \sqrt{r})} \sin(1 / r)\] $\qquad$ where $\; r = \sqrt{x^2 + y^2} $

    と表される。(8.4.c)から指数部分以降を省くと線形自励形となり [1] - 6 で見たように原点から放射状に伸びる半直線である。 では、(8.4.c)自体はどのようになるだろうか?

    (8.4.c)のベクトル場を示すと以下のようになる。

    ベクトル場のグラフは各線分の向きを矢印ではなく、赤色のドットで代用している。 ベクトルは放射状であるが、$\alpha(r) = 0$ となる円周 $r = 1 / k\pi$ を境に向きが変わるのが確認できると思う。 一番外側の黄色の円周を例に取ると、円周の外側の各点(Mとする)は原点から放射状に遠ざかるベクトルである。 黄色の円周と1つ内側の青色の円周までの各点(Nとする)は、放射状に原点に向かうベクトルである。

    話を方程式(8.4.a)に戻す。 方程式(8.4.a)において、黄色の円周より外側の点Mのベクトルは、 原点から点Mまでの距離を半径とする円に沿ったベクトルと原点から放射状に遠ざかるベクトルの和である。 解軌道は、円周に沿いながら離れていく軌道になる。
    一方、黄色の円周と1つ内側の青色の円周間の点Nのベクトルは、 原点から点Nまでの距離を半径とする円に沿ったベクトルと放射状に原点に向かうベクトルの和である。 解軌道は、青色の円周に沿って向かいながら、やがて青色の円周に巻き付いて行く軌道になる。
    更に、青色の円周からその内側の黄色の円周までの各点は、 原点から点までの距離を半径とする円に沿ったベクトルと放射状に外側に向かうベクトルの和になるので、 解軌道は、青色の円周に沿って向かいながら、やがて青色の円周に巻き付いて行く軌道になる。

    原点に近づくに連れて、円周 $1 / k\pi$ は無限に出現するが、 上記の説明の繰返しになるため、$k$が奇数の場合、極限円は不安定となり、 偶数の場合は安定となる。

  9. Hamilton型の例 1  Cassini(カシニ)の燈形

    例)

    \[ \frac{dx}{dt} = y(1 + x^2 + y^2)\] \[ \frac{dy}{dt} = x(1 - x^2 - y^2)\]

    ハミルトン型の方程式は以下の特徴を持つ。

    \[ \frac{dx}{dt} = F_y(x, y)\] \[ \frac{dy}{dt} = -F_x(x, y)\] \[ 解析解は、以下の陰関数となる \] \[ F(x, y) = C\]

    解析解が陰関数となるので、直接、一般解を求積しても、描画は以外とやっかいである。 媒介変数を使って表示するという手もありだが、関数に手を加えず描きたいと横着な考え方をすると、 描画エリアを格子状に分割して、各格子の4辺と関数が交わるか調べながら 近似的にグラフをプロットしていくというのが、簡単な方法ではないかと思う。

    微分方程式の解軌道は、解析解の曲線族をプロットすることになるので、 Hamilton型の場合は、ダイレクトに陰関数のグラフを描画することになる。

    Cassiniの燈形の解析解は以下の通りである。

    \[ F(x, y) = \frac{1}{4}(1 + x^2 + y^2)^2 - x^2 = C\]

    平衡点は原点と、$(\pm1, 0)$の3点である。原点は鞍状点になり、$(\pm1, 0)$は渦心点となる。

    実際にプロットしてみると、解軌道は全て周期解である。 陰関数のグラフ Cassiniの燈形 と一致することが確認できる。

    陰関数のグラフの定数が0.25のとき、レムニスケートになる。 [4] - 6 でレムニスケートを紹介したが、一般解は極限軌道であったことが大きく異なる。

  10. Hamilton型の例 2  Descartes(デカルト)の葉線

    例)

    \[ \frac{dx}{dt} = 2x - 3y^2\] \[ \frac{dy}{dt} = -2y + 3x^2\]

    解析解は以下の通りである。

    \[ F(x, y) = -x^3 + 2xy - y^3 = C\]

    平衡点は原点と、$(2/3, 2/3)$の2点である。原点は鞍状点になり、$(2/3, 2/3)$は渦心点となる。

    陰関数のグラフ で定数が0の場合、デカルトの正葉線がプロットされる。

  11. Hamilton型の例 3  Descartes(デカルト)の葉線の派生形

    例)

    \[ \frac{dx}{dt} = 2x^2 - 5y^4\] \[ \frac{dy}{dt} = -4xy + 5x^4\]

    解析解は以下の通りである。

    \[ F(x, y) = -x^5 + 2x^2y - y^5 = C\]

    平衡点は原点と、第1, 3象限(それぞれ点対称の関係)の3点である。原点は鞍状点になり、1, 3象限は渦心点となる。

    陰関数のグラフ で定数が0の場合、デカルトの正葉線を2つ組み合わせたような曲線がプロットされる。

    確認のために、曲面$F(x, y)$のグラフを描くと以下のようになる。 陰関数のグラフが曲面の等高線になっていることが分かる。

  12. Hamilton型の例 4  三角関数の周期性が反映される例

    例)

    \[ \frac{dx}{dt} = \sin(y)\] \[ \frac{dy}{dt} = \sin(x)\]

    解析解は以下の通りである。

    \[ F(x, y) = \cos(y) - \cos(x) = C\]

    平衡点は $(p\pi, q\pi)$ である($p,q$ は整数)。 $p + q$ が偶数のときは鞍状点になり、奇数のときは渦心点になる。

    陰関数のグラフ からは、$p + q$ が奇数で、$p$ が奇数のときは凸の等高線になり、 $p$が偶数のときは凹の等高線になることを示している。

    曲面$F(x, y)$のグラフを描くと以下のようになる。 これ、台所で見かける何かに似ていると思いませんか?

[99] ユーザー定義

[描画へ]をクリックすると、以下の画面が表示されます。 任意の方程式を入力して、その解軌道を確認することができます。 STARTボタンを押すと入力ボックスの下に入力した式がTeX形式で表示され、描画が開始されます。 入力された式の評価には Math.js を利用しています。 エラーメッセージが表示された場合は、入力した式を確認して、適宜修正して下さい。
入力できる方程式は自励形に限ります(非自励形については、現在確認中)。

以下に、入力例を示します。

  1. 2次元1次線形自励系の例です。

  2. 連立方程式にする必要がない場合の例です。$dx/dt$ を1とします。

  3. べき乗は様々な入力が可能ですが、単純なものであれば、x^2で十分です。勿論、x*xでもよいです。 2乗くらいならx*xで問題ないですが、4乗、5乗となってくるとx^4などの方が入力が楽です。
    pow()も使えます。 pow(x^2 + y^2, 2)などはTeXで表示したとき、$(x^2 + y^2)^2$と表示されるので違和感がありませんが、 pow(x, 5)と入力すると、表示は$(y)^5$となります。計算自体は正しく実行されるので問題はないのですが、 違和感を感じられる場合は、単純にy^5と入力して下さい。

  4. 式が長くなる場合は、以下のように折り返して入力できます(ただし、異常に長い場合は画面のレイアウトが崩れます。 崩れても計算に影響はありませんが、非常に見づらくなります)。 また、式の一部が繰り返されるときは、where, $r = $ の入力ボックスに式を入力して、 方程式で繰り返される部分をrで置換することで、入力の手間を省くことができます。
    必要ない場合は、1.~3.のように何も入力する必要はありません。

  5. 入力された数式の解析はmath.jsを利用しているため、 pow()も含めて、sin(), cos()などの数学関数はそのまま入力できます。
    逆に、入力ボックスにmath.pow(), math.sin(), math.cos()などを含めた数式が入力された場合は、 エラーになり、描画は実行されません。
    指数については、e^xは$e^x$として評価されます。 exp(x)と入力した場合も、TeXの表示は$\exp(x)$になりますが、同じように評価されます。 指数部が複雑な場合はe^(x + 1)或いはexp(x + 1)のように括弧をつけて下さい。

== 注釈 ==

[1.]

このページでは、Local Storageを利用してセレクトボックスで選択されたリストの位置を保持する機能を持ちます。 選択項目の数が相当な量のため、ブラウザの更新ボタンが押された場合や、ページを閉じた後に再度ページを開いたとき、 いちいち先頭のリストから目的のリストまで辿るのは大変です。

そこで、以下に示すように選択されたリストの位置を保持するためのチェックボックスを設置しました。

チェックを入れてから特定の項目を選択すると、 Local StorageにチェックボックスのON/OFFと選択されたリストの番号が保存されます。 ブラウザの更新ボタンの押下やページを閉じてもデータが保存されるため、 ページが表示されると同時に目的の項目が選択された状態になります。

ただし、チェックを外すとLocal Storageに保存されていたデータは、keyと共に削除されるため、 ブラウザの更新ボタンが押された場合や、ページを閉じた後に再度ページを開いたとき、リストの先頭が選択された状態に戻ります。

Local Storageに保存されるデータ ( key, vakue ) は次の2点です。

keyvalue
kzDeq_StorageFlg"true" or "false"
kzDeq_SelectedType"09&09" などの5文字の文字列

Local Storageに保存されているデータは、お使いのブラウザの設定画面から確認することができます。 戻る