自習室

こもります

JavaScriptでニュートン法を用いて三次方程式を解いて、ついでにグラフも描く

初めに

完成品

ここで動きます

動機

そもそもの目的は、楕円の法線を求めるために四次方程式を解くことだったのですが、まずは手始めに三次方程式を解くことにしました。結果、高校時代に勉強した内容を思い出して、懐かしいような気分に浸れて、変な楽しさがありました(笑)

解の公式を使えば良いのでは?

三次方程式も四次方程式も解の公式が知られていますが、たとえば三次方程式を解く「カルダノの方法」だと、

{ \displaystyle
B = x^{3} - 15x - 4
}

という式を解いて x = 4 を算出することが出来ないなど、汎用性の乏しさが知られています。

三次方程式 - Wikipedia

ニュートン法を使う

この手の計算では一般的に使われているらしいニュートン法を利用する事にしました。

Newton 法による方程式の近似解法

三次方程式は実数解を最低ひとつは持つことが明かですので、まずはニュートン法で実数解をひとつ求めて、残りの解(虚数解も含む)を求めるような関数を作ります。

コード

コードはこちらに上げてあります。

izmhr/newtonCubicEquation · GitHub

場合分け

ポイントだけ解説します

計算をやりやすくするために、特徴的な係数や特徴的な解が得られる場合をあらかじめ場合分けします。やり方にはいろいろあると思いますが、今回のプログラムで採用している場合分けについて紹介します。

{
f(x) = ax^{3} + bx^{2} + cx + d = 0
}

を解きます。

解に0 が含まれる場合

言い換えると、定数項 {d=0}の場合。

{ \displaystyle
f(x) = ax(x^{2} + \frac{b}{a}x + \frac{c}{a}) = 0
}

を解くことになります。 これは一般的な二次方程式の判別式を用いて、残りの2解がどうなるかを場合分けします

{
D = b^{2} - 4ac
}

とすると、

  • { D \gt 0} なら、0に加えて異なる実数解が2つ
  • { D = 0} なら、0に加えて、実数の重解1組
  • { D \lt 0} なら、0に加えて、共役虚数解1組

がそれぞれ存在します。

解に0が含まれない場合は、ニュートン法を用いて、1つ目の実数解を求める

が、そこであたえた近似開始値が解そのものだった場合

ニュートン法の近似開始値 {x_k} に対して

{
ax_k^{3} + bx_k^{2} + cx_k + d = 0
}

この{x_k} をそのまま解の1つ目の実数解として利用する

近似開始値が、極値をとるxであった場合

この場合、ニュートン法が収束しなくなるので、別の近似開始値 {x_k} を使うよううながす。私のプログラムでは0.1加えるようにしているが、万能ではないのでご注意です。

ニュートン法で1つ目の実数解が求まった後

ここでも判別式を利用する。ここで利用する解のパタンを判別する式については、こちらのサイトが詳しいです

http://www2.odn.ne.jp/~aai55890/donnwa2/sanzihannbetu.htm

 {
D_1 = b^{2} - 3ac
}

 {
D_2 = 27a^{2}b^{2} - 18abcd - b^{2}c^{2} + 4b^{3}d + 4ac^{3}
}

  • { D_1} は、{ f^{'}(x)} の判別式。グラフの形状の把握に使う。
  • { D_2} は、{f(x)}極値 {f(\alpha), f(\beta)} に対して、{f(\alpha)f(\beta)} で算出される値で、これも判別式として利用する。x軸との交差回数の判定のために使います

詳細は先ほどのサイトに。

{ D_1 \gt 0} かつ { D_2 \lt 0} のとき

ことなる実数解3つが存在する。1つ目の実数解{x_0}を用いると、解きたい方程式は

{
f(x) = ax^{3} + bx^{2} + cx + d = ax(x - x_0)(x^{2} + Ax + B) = 0
}

と書き直すことが出来、二次方程式を解けば良いことになる。

{ D_1 \gt 0} かつ { D_2 = 0} のとき

このとき、異なる実数解2つで片方は重解になる。 ここで、先のニュートン法で求まった実数解が重解の一部である場合と、独立した解である場合の2通りがあり得る。

ここで、{ x^{2} + Ax + B = 0} が、重解を持てば、{x_0} は単独の解であったということになりますが、{x_0}ニュートン法で求めた解であるため、導き出されるAやBも、誤差を含む値になります。従って、 { x^{2} + Ax + B = 0} が重解を持つかどうかの判別は、

{ A^{2} - 4B = 0}

ではなく

{|A^{2} - 4B| \lt Thresh}

(ここで Threshは、十分に小さい値) という形で行うことになります。 これが満たされれば、{ x^{2} + Ax + B = 0} は重解を持ち、元の{x_0} に加えてもう1解算出すれば良いです。

満たされていない場合は、

{
f(x) = ax^{3} + bx^{2} + cx + d = ax(x - x_0)^{2}(x - x_1) = 0
}
{ x_1 = -\frac{d}{x_0^{2}a}}

が、もう一つの解になります。

それ以外の場合

実数解は初めの1つだけで、残りは共役な虚数解を持つことになります

ウェブアプリとして仕上げる

機能

機能として、以下の様な物を実装しています

  • 主要機能
  • 可視化関連の機能
    • 与えられた三次関数のグラフを描画
    • 実数解をプロット
    • グラフのズームの調整
  • (おまけ) 手っ取り早く試せるよう、サンプルの関数を5つプリセット

CreateJS を利用

グラフの描画は、Canvasに対してCreajteJSを使って行っています。D3などの利用も考えましたが、最近よく使っていて手が慣れてたのでCreateJSを利用しました。

さいごに

  • 算術
  • Canvasの描画
  • ツールとして機能を設計してまとめ上げる
  • htmlのインプットタグなどを使いこなす

など、結構やることが多くて、勉強になりました。インプットタグを使うのが意外としんどかったです。

実際は虚数解が必要になるようなケースはあまりないので、正しい制約下で必要な実数解を最短で算出することが求められるようなことの方が多いでしょう。まだほとんど知らないのですが、次は非線形最適化のツールなど使ってみようと思っています。

 {
Y(n) = a_nX(n) + a_{n-1}X(n-1)
}