JavaScriptで連立方程式を解きたかった話 #JavaScript

経緯

弊社のサービスでmockmockというIoTに関わるエンジニア向けのサービスがあるのですが、その中で「グラフをグニグニ動かせるUIが欲しい」という要望を振られました。

確かにAdobe AfterEffectだったりMayaだったり、そういうツールにはグラフをグニグニ動かすUIっていうのは欠かせないものですが、我流で作ってみました。

点Aと点Bの２点を結ぶとき、以下の条件をみたすようにしました。

点Aを通過する

点Bを通過する

点Aの接線の傾きを指定する

点Bの接線の傾きを指定する

そのような三次関数で点Aと点Bの間を補完する。

方針

三次関数は一般式で

$y = ax^3 + bx^2 + cx + d$

で表すことができます。今回の目標は$(a,b,c,d)$の組を見つけること。

ひとまず条件を数式化。

$y$を$x$で微分すると、

$\frac{dy}{dx} = 3ax^2 + 2bx + c$

点A、点Bをそれぞれ、$(x_1, y_1)$,$(x_2, y_2)$、点A、点Bでの傾きをそれぞれ$\alpha$,$\beta$とおくと、以下の通り。

$x_1^3a + x_1^2b + x_1c + d = y_1$
$x_2^3a + x_2^2b + x_2c + d = y_2$
$3x_1^2a + 2x_1b + c = \alpha$
$3x_2^2a + 2x_2b + c = \beta$

この連立方程式を$(a,b,c,d)$について解こうとしたわけです。

解く

最初は手書きでといて、その解をコードに落とせばいいや、と思ってたのですが、メンドイ。計算なんてものはコンピュータにやらせりゃいいのです。計算機としての本分でしょう。

方針としてクラーメルの公式を使いました。

ひとまず、行列化。

\begin{equation}
\begin{pmatrix}
x_1^3 & x_1^2 & x_1 & 1 \\
x_2^3 & x_2^2 & x_2 & 1 \\
3x_1^2 & 2x_1 & 1 & 0 \\
3x_2^2 & 2x_2 & 2 & 0
\end{pmatrix}
\begin{pmatrix}
a \\
b \\
c \\
d
\end{pmatrix}
= 
\begin{pmatrix}
y_1 \\
y_2 \\
\alpha \\
\beta
\end{pmatrix}
\end{equation}

クラーメルの公式を使うと、$A$の行列式を$|A|$として

\begin{equation}
a = \frac{
\begin{vmatrix}
y_1 & x_1^2 & x_1 & 1 \\
y_2 & x_2^2 & x_2 & 1 \\
\alpha & 2x_1 & 1 & 0 \\
\beta & 2x_2 & 2 & 0
\end{vmatrix}
}
{
\begin{vmatrix}
x_1^3 & x_1^2 & x_1 & 1 \\
x_2^3 & x_2^2 & x_2 & 1 \\
3x_1^2 & 2x_1 & 1 & 0 \\
3x_2^2 & 2x_2 & 2 & 0
\end{vmatrix}
},

b = \frac{
\begin{vmatrix}
x_1^3 & y_1 & x_1 & 1 \\
x_2^3 & y_2 & x_2 & 1 \\
3x_1^2 & \alpha & 1 & 0 \\
3x_2^2 & \beta & 2 & 0
\end{vmatrix}
}
{
\begin{vmatrix}
x_1^3 & x_1^2 & x_1 & 1 \\
x_2^3 & x_2^2 & x_2 & 1 \\
3x_1^2 & 2x_1 & 1 & 0 \\
3x_2^2 & 2x_2 & 2 & 0
\end{vmatrix}
},

c = \frac{
\begin{vmatrix}
x_1^3 & x_1^2 & y_1 & 1 \\
x_2^3 & x_2^2 & y_2 & 1 \\
3x_1^2 & 2x_1 & \alpha & 0 \\
3x_2^2 & 2x_2 & \beta & 0
\end{vmatrix}
}
{
\begin{vmatrix}
x_1^3 & x_1^2 & x_1 & 1 \\
x_2^3 & x_2^2 & x_2 & 1 \\
3x_1^2 & 2x_1 & 1 & 0 \\
3x_2^2 & 2x_2 & 2 & 0
\end{vmatrix}
},

d = \frac{
\begin{vmatrix}
x_1^3 & x_1^2 & x_1 & y_1 \\
x_2^3 & x_2^2 & x_2 & y_2 \\
3x_1^2 & 2x_1 & 1 & \alpha \\
3x_2^2 & 2x_2 & 2 & \beta
\end{vmatrix}
}
{
\begin{vmatrix}
x_1^3 & x_1^2 & x_1 & 1 \\
x_2^3 & x_2^2 & x_2 & 1 \\
3x_1^2 & 2x_1 & 1 & 0 \\
3x_2^2 & 2x_2 & 2 & 0
\end{vmatrix}
}
\end{equation}

コードにするべ

JavaScriptで行列計算をやろうと思ったら、まー、math.jsでしょう。
math.jsがある状態でやってみた。
JavaScriptといいながら、CoffeeScriptで。

x1 = 0
y1 = 0
alpha = 0
x2 = 100
y2 = 100
beta = 0

# Av = bの形の連立方程式を解く
# @param num [Integer] 変数の数
# @param coefficients [Array<Array<Float> >] 係数行列A: num x num 行列
# @param results [Array<Float>] 右辺のベクトルb: num x 1 行列
solveSimultaneousEquation = (num, coefficients, results) ->
  # 係数行列の行列式の計算
  coefficientsMat = math.matrix(coefficients)
  coefficientsDet = math.det(coefficientsMat)

  solutions = []
  for j in [0..num-1]
    replacedCoefficients = coefficients.concat()
    for k in [0..num-1]
      replacedCoefficientsRow = replacedCoefficients[k].concat()
      replacedCoefficientsRow[j] = results[k]
      replacedCoefficients[k] = replacedCoefficientsRow
    replacedCoefficientsDet = math.det(math.matrix(replacedCoefficients))
    solutions.push replacedCoefficientsDet / coefficientsDet

  return solutions

coefficients = [
  [ 1,x1, Math.pow(x1, 2),     Math.pow(x1, 3) ]
  [ 1,x2, Math.pow(x2, 2),     Math.pow(x2, 3) ]
  [ 0, 1,          2 * x1, 3 * Math.pow(x1, 2) ]
  [ 0, 1,          2 * x2, 3 * Math.pow(x2, 2) ]
]
results = [
  y1
  y2
  alpha
  beta
]

solutions = solveSimultaneousEquation(4, coefficients, results)

けっこうキモなのは列の順番を逆にしているところ。

列の順番をそのまま[x1^3, x1^2, x1, 0]みたいにしてしまうと、$x_1$とかが大きい値になった場合、$x_1^3$がデカくなって、行列式の計算をする際、情報落ちで誤差がひどくて、全然違う結果になってしまいます。情報落ちの誤差の影響が顕在化したのは初めての経験でした。

まとめ

やっぱり数学楽しい。have a nice Math!¹

私は黒い三角定規の一員ではございません ↩