Posted at
FusicDay 14

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

More than 1 year has passed since last update.


経緯

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

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

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



  1. 点Aを通過する

  2. 点Bを通過する

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

  4. 点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$とおくと、以下の通り。


  1. $x_1^3a + x_1^2b + x_1c + d = y_1$

  2. $x_2^3a + x_2^2b + x_2c + d = y_2$

  3. $3x_1^2a + 2x_1b + c = \alpha$

  4. $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!1





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