算数 数学

音源定位　update 20180429

音波を発する原因を音源と呼びます。
例えば、話し声は音波で、しゃべっている本人を音源ということができます。また、工事現場では騒音は音波で、その騒音を発している重機も音源といえます。

そんな音源に対して、複数のマイクロホン（以下マイクと呼びます）を並べて、音波が各マイクに到達する時間差から、音源の位置を推定することができます。これを音源定位といいます。

初めから立体的な配置を考えると、とても複雑になるので、音源もマイクも平面上にあると考えましょう。立体的な配置を考えるのは平面の話が済んでからでも遅くはないと思います。

マイクが1個だと、音源の有る・無しと音の強弱がわかります。ですが、どちらの方向から音がやってくるのかがわかりません。

マイク２個だと、人間の耳を連想するので、なんとなく音のやってくる方向までは分かるんじゃないのかと思います。

１個の音源とマイク１、マイク２が平面上に配置されている場合で計算してみましょう。

音源１個の対してマイクが２個の場合の関係式を求める

（x0 ，y0 ）：音源の座標
（x1，y1）：マイク1の座標
（x2，y2）：マイク2の座標
v：音速
L1：音源からマイク1までの距離，
L2：音源からマイク2までの距離
t1：マイク1の音波到達時刻，
t2：マイク2の音波到達時刻

t2＝L2/v　（t2：マイク2の音波到達時刻，L2：音源からマイク2までの距離，v：音速）
t1＝L1/v　（t1：マイク1の音波到達時刻，L1：音源からマイク1までの距離，v：音速）
L2 ＝SQRT((x2-x0)^2+(y2-y0)^2) 　　（（x0，y0）：音源の座標）
L1 ＝SQRT((x1-x0)^2+(y1-y0)^2) 　　（（x2，y2）：マイク2の座標、 （x1，y1）：マイク1の座標）
Δ＝(t2-t1)＝(L2-L1)/v

(t2-t1)*v＝L2-L1＝SQRT((x2-x0)^2+(y2 -y0)^2) -SQRT((x1-x0)^2+(y1-y0)^2)

(t2-t1)^2*v^2＝(x2-x0)^2+(y2 -y0)^2+(x1-x0)^2+(y1-y0)^2+2*SQRT((x2-x0)^2+(y2-y0)^2) *SQRT((x1-x0)^2+(y1-y0)^2)

(t2-t1)^2*v^2−((x2−x0)^2 ＋ (y2 −y0)^2 ＋(x1−x0)^2 ＋ (y1 −y0)^2)＝2*SQRT((x2−x0)^2 ＋ (y2 −y0)^2) *SQRT((x1−x0)^2 ＋ (y1 −y0)^2)

左辺^2＝ ((t2-t1)^2*v^2-((x2-x0)^2+(y2-y0)^2+(x1-x0)^2+(y1-y0)^2))^2

右辺^2＝ 4*((x2-x0)^2+(y2-y0)^2)*((x1-x0)^2+(y1-y0)^2

したがって
0=左辺^2-右辺^2
=((t2-t1)^2*v^2-((x2-x0)^2+(y2-y0)^2+(x1-x0)^2+(y1-y0)^2))^2-4*((x2-x0)^2+(y2-y0)^2)*((x1-x0)^2+(y1-y0)^2)

∴
+4*((x2-x1)^2-(t2-t1)^2*v^2)*x0^2
+4*(-(x2-x1)*((y2-y1)*(y2+y1)+(x2-x1)*(x2+x1))+(t2-t1)^2*(x2+x1)*v^2)*x0
+8*(x2-x1)*(y2-y1)*x0*y0
+4*((y2-y1)^2-(t2-t1)^2*v^2)*y0^2
+4*(-(y2-y1)*((y2-y1)*(y2+y1)+(x2-x1)*(x2+x1))+(t2-t1)^2*(y2+y1)*v^2)*y0
+((y2-y1)*(y2+y1)+(x2-x1)*(x2+x1))^2-2*(t2-t1)^2*v^2*(y2^2+y1^2+x2^2+x1^2)+(t2-t1)^4*v^4
=0

音源座標は、マイク1(x1，y1)と、マイク2(x2，y2)と音波到達時刻の差(t2−t1)、および音速vによる係数を持つ2次曲線上にある。（この場合の2次曲線は双曲線）

ページのトップへ

----- ここまで　20180422 -----

音源位置(x0, y0)を推定する計算式を求める。

各マイクの位置を(x1, y1), (x2, y2), (x3, y3)、音源位置を(x0, y0)、
音の発生時刻をt0（未知）、各マイクに音が到達する時刻を t1, t2, t3 とする。
・マイク１，２に音が到達する時間差(t2-t1)から次式が成り立つ。

SQRT((x2-x0)^2+(y2-y0)^2) −SQRT((x1-x0)^2+(y1-y0)^2) = v*(t2-t1)　…�@

v : 音速（既知、気温により決まる） 　

同様に、マイク２，３に音が到達する時間差(t3-t2)から次式が成り立つ。

SQRT((x3-x0)^2+(y3-y0)^2) −SQRT((x2-x0)^2+(y2-y0)^2) = v*(t3-t2)　…�A

しかし、ここでは音源位置を求めるために以下の手順を採用する。
・音源から各マイクまでの距離と音の到達時間より次式が成り立つ。

SQRT((x1-x0)^2+(y1-y0)^2) = v*(t1-t0) 　…�C

SQRT((x2-x0)^2+(y2-y0)^2) = v*(t2-t0) 　…�D

SQRT((x3-x0)^2+(y3-y0)^2) = v*(t3-t0) 　…�E

・�C〜�Eを2乗すると、
(x1-x0)^2 + (y1-y0)^2 = (v*(t1-t0))^2　…�F

(x2-x0)^2 + (y2-y0)^2 = (v*(t2-t0))^2 　…�G

(x3-x0)^2 + (y3-y0)^2 = (v*(t3-t0))^2 　…�H

・�Fから�Gを引くと、
(x1-x0)^2 + (y1-y0)^2 - (x2-x0)^2 - (y2-y0)^2 = (v*(t1-t0))^2 - (v*(t2-t0))^2

展開整理して、
(x1^2+y1^2-x2^2-y2^2)-2*x0*(x1-x 2) -2*y0*(y1-y2) = (v^2)*(t1^2-t2^2)-2*(v^2)*t0*(t1-t2)　…�I

・同様に、�Fから�Hを引き、展開整理すると、
(x1^2+y1^2-x3^2-y3^2)-2*x0*(x1-x 3) -2*y0*(y1-y3) = (v^2)*(t1^2-t3^2)-2*(v^2)*t0*(t1-t3) 　…�J

・�I、�Jを変形して、
-2*(y1-y2)*y0+2*(v^2)*(t1-t2)*t0 = (v^2)*(t1^2-t2^2)-(x1^2+y1^2-x2^2-y2^2)+2*(x1-x2)*x0 　…�K

-2*(y1-y3)*y0+2*(v^2)*(t1-t3)*t0 = (v^2)*(t1^2-t3^2)-(x1^2+y1^2-x3^2-y3^2)+2*(x1-x3)*x0　…�L

・ここで、

　　a1=-2*(y1-y2),
　　b1= 2*((v^2)*(t1-t2),
　　c1= (v^2)*(t1^2-t2^2)-(x1^2+y1^2-x2^2-y2^2),
　　d1= 2*(x1-x2)

　　a2= -2*(y1-y3),
　　b2= 2*(v^2)*(t1-t3),
　　c2= (v^2)*(t1^2-t3^2)-(x1^2+y1^2-x3^2-y3^2),
　　d2= 2*(x1-x3)

とおくと、�K、�Lは
　　a1*y0 + b1*t0 = c1 + d1*x0　…�M
　　a2*y0 + b2*t0 = c2 + d2*x0　…�N となる。

・これをy0, t0に関する連立方程式と考えて解くと、
　　y0 = A*x0 + B　…�O
　　t0 = C*x0 + D　…�P

ここで、A = (b2*d1-b1*d2)／(a1*b2-a2*b1),
B = (b2*c1-b1*c2)／(a1*b2-a2*b1)
　　　 C = (a1*d2-a2*d1)／(a1*b2-a2*b1),
D = (a1*c2-a2*c1)／(a1*b2-a2*b1)

・この y0、t0を�Fに代入すると、x0に関する次の2次方程式を得る。

　　(x1-x0)^2 + (y1-(A*x0+B))^ 2 = (v*(t1-(C*x0+D)))^2

x0について整理すると、

　　E*x0^2 + F*x0 + G = 0

ここで、
　　　E = 1+A^2-(v^2)*(C^2)
　　　F = -2*x1-2*(y1-B)*A +2*(v^2)*(t1-D)*C
　　　G = x1^2+(y1-B)^2-(v^2)*(t1-D)^2

・この2次方程式を解くと音の発生位置x0が得られ、
更に�Oよりy0が求められる。

　　x0 = (-F±SQRT(F^2 - 4*E*G))／(2*E)　　（但し、F^2 - 4*E*G ≧ 0 のとき）

（注）x0に関する2次方程式からは通常解が２個得られるが、
　　　このx0を�Pに代入して 得られるt0について、
　　　t0 ＜ t1, t2, t3 を満たさないものは除外する。

ここでx0(+)＝ (-F+SQRT(F^2 - 4*E*G))／(2*E)
　　　x0(−)＝ (-F−SQRT(F^2 - 4*E*G))／(2*E)
とすると
　　　　　y0(+)＝A*x0(+) + B
　　　　　y0(−)＝A*x0(−) + B
となる。

音源は１個ですが、計算式からは2個求まることになります。
どちらか一方は虚像ということです。

ページのトップへ

----- ここまで　20180423 -----

マイクを４個、５個と増やしていっても、これまでやってきた連立方程式を解く方法では、音源のx座標を求める式の次数を２次から１次に下げることはできないようです。
２次のままだと、常に虚像が現われてしまいます。

そこで、マイクを６個準備して、３個づつ２個のグループでそれぞれ音源の位置を計算するとします。
今、マイクＭ１、Ｍ２、Ｍ３のグループで求められる音源のX座標を、x=b/a（真の音源のx座標）とx＝d/c（虚像）、
マイクＭ４、Ｍ５、Ｍ６のグループで求められる音源のX座標を、x=b/a（真の音源のx座標）とx=-f/e（虚像）とします。

３個のマイクで求められる音源座標は、常に、真の座標と虚像の座標ののペアになります。どのマイクのグループを作っても、それから求められる２個の座標のうち１つは真の音源座標ですから、同じ解を１個持つ２次の連立方程式を考えることができます。

そこで、
1個の共通の根x=b/aを持つ2つの2次方程式を考える。
(ax-b)(cx-d)=0
(ax-b)(ex-f)=0
展開すると
acx＾2-adx-bcx+bd=0
aex＾2-afx-bex+bf=0
a≠0、c≠0、e≠0のとき
x＾2の項を消去する
a＾2cex＾2-a＾2dex-abcex+abde=0・・・�M
a＾2cex＾2-a＾2cfx-abcex+abcf=0・・・�N
�M-�Nは以下のとおり
a＾2(cf-de)x-ab(cf-de)=0
cf≠deのとき、（これは2つの2次方程式は異なることを意味する）
a＾2x-ab=0
a≠0だから
ax=b
∴　x=b/a　→　真の音源のX座標が導かれました。

マイクを６個使うと、虚像を消せそうなことがわかりました。
しかし、マイクを６個使っても、一直線に並べてしまうと、虚像が同じ場所にできる（x座標が同じになる）ので意味がなくなります。
虚像が異なる場所にできるようにマイクの配置を考えなくてはなりません。

一直線に並べなければ、1個のマイクを共通に使うことで５個で良いことになります。

この場合、Ｍ１、Ｍ２、Ｍ３でひとつのグループ、Ｍ３、Ｍ４、Ｍ５でもうひとつのグループとして計算をします。
ひとつのグループ内のマイクの並び方は一直線です。

３個のマイクの並び方を直線ではなく不規則にすると、４個で良いことになります。

ページのトップへ

----- ここまで　20180426 -----

マイク４個の位置関係を以下のように設定します。

それぞれのマイクのx-y座標を次のようにします。この時の音速を、334m/ｓとして、マイクでの測定精度を1000分の1秒とします。

Ｍ１：（ｘ，ｙ）＝（0m，1m）、測定時刻＝12.700秒
Ｍ２：（ｘ，ｙ）＝（3000m，-1m）、測定時刻＝8.985秒
Ｍ３：（ｘ，ｙ）＝（1000m，30m）、測定時刻＝10.720秒
Ｍ４：（ｘ，ｙ）＝（2000m，-30m）、測定時刻＝9.5530秒

Ｍ１〜Ｍ２〜Ｍ３の組み合わせによる音源座標の計算は

及び

となります。
この段階では、どちらが真の音源座標かわかりません。

次にＭ１〜Ｍ２〜Ｍ４の組み合わせによる音源座標を計算します。

及び

このように計算結果を並べてみると、真の音源座標はなんとなく（ｘ，ｙ）＝（3000ｍ、3000ｍ）と予想できます。

実際、音源座標は

に設定していました。

ここで、虚像を消して一気に真の音源座標を求めてみましょう。

マイクの組み合わせＭ１〜Ｍ２〜Ｍ３による計算の地中で使用した係数　Ｅ１，Ｆ１，Ｇ１と、
Ｍ１〜Ｍ２〜Ｍ４での係数　Ｅ２，Ｆ２、Ｇ２を用いて

ｙ座標を２個計算しているのは、それぞれのマイクの組み合わせの途中計算で使用した係数Ａ，Ｂを使用した結果で、どちらを使用しても良いし、平均値をとっても良いのではと思います。

ページのトップへ

----- ここまで　20180428 -----

「音源１個にマイク２個の場合」の章と、「音源１個のマイク３個の場合」の章で使用した計算式の形が違うことに気が付きましたか？

「音源１個にマイク２個の場合」の章で使った計算式は

「音源１個のマイク３個の場合」の章で使った計算式は

です。
一見、音源が音波を発した時刻t0を含まない上の式の方が良いように思えます。

そこで、こちらの方で計算してみましょう。
こちらの式からは、
+4*((x2-x1)^2-(t2-t1)^2*v^2)*x0^2
+4*(-(x2-x1)*((y2-y1)*(y2+y1)+(x2-x1)*(x2+x1))+(t2-t1)^2*(x2+x1)*v^2)*x0
+8*(x2-x1)*(y2-y1)*x0*y0
+4*((y2-y1)^2-(t2-t1)^2*v^2)*y0^2
+4*(-(y2-y1)*((y2-y1)*(y2+y1)+(x2-x1)*(x2+x1))+(t2-t1)^2*(y2+y1)*v^2)*y0
+((y2-y1)*(y2+y1)+(x2-x1)*(x2+x1))^2-2*(t2-t1)^2*v^2*(y2^2+y1^2+x2^2+x1^2)+(t2-t1)^4*v^4
=0
が導かれました。

これは、マイクＭ１〜Ｍ２を焦点とする双曲線を表していますが、
同様にＭ２〜Ｍ３、Ｍ３〜Ｍ１でも双曲線が描けます。

まとめると、

ただしA〜Sの係数は以下の通り、

A＝+4*((x2-x1)^2-(t2-t1)^2*v^2)
B＝+4*(-(x2-x1)*((y2-y1)*(y2+y1)+(x2-x1)*(x2+x1))+(t2-t1)^2*(x2+x1)*v^2)
C＝+8*(x2-x1)*(y2-y1)
D＝+4*((y2-y1)^2-(t2-t1)^2*v^2)
E＝+4*(-(y2-y1)*((y2-y1)*(y2+y1)+(x2-x1)*(x2+x1))+(t2-t1)^2*(y2+y1)*v^2)
F＝+((y2-y1)*(y2+y1)+(x2-x1)*(x2+x1))^2-2*(t2-t1)^2*v^2*(y2^2+y1^2+x2^2+x1^2)+(t2-t1)^4*v^4

G＝+4*((x3-x1)^2-(t3-t1)^2*v^2)
H＝+4*(-(x3-x1)*((y3-y1)*(y3+y1)+(x3-x1)*(x3+x1))+(t3-t1)^2*(x3+x1)*v^2)
I＝+8*(x3-x1)*(y3-y1)
J＝+4*((y3-y1)^2-(t3-t1)^2*v^2)
K＝+4*(-(y3-y1)*((y3-y1)*(y3+y1)+(x3-x1)*(x3+x1))+(t3-t1)^2*(y3+y1)*v^2)
L＝+((y3-y1)*(y3+y1)+(x3-x1)*(x3+x1))^2-2*(t3-t1)^2*v^2*(y3^2+y1^2+x3^2+x1^2)+(t3-t1)^4*v^4

M＝+4*((x3-x2)^2-(t3-t2)^2*v^2)
N＝+4*(-(x3-x2)*((y3-y2)*(y3+y2)+(x3-x2)*(x3+x2))+(t3-t2)^2*(x3+x2)*v^2)
O＝+8*(x3-x2)*(y3-y2)
P＝+4*((y3-y2)^2-(t3-t2)^2*v^2)
R＝ +4*(-(y3-y2)*((y3-y2)*(y3+y2)+(x3-x2)*(x3+x2))+(t3-t2)^2*(y3+y2)*v^2)
S＝+((y3-y2)*(y3+y2)+(x3-x2)*(x3+x2))^2-2*(t3-t2)^2*v^2*(y3^2+y2^2+x3^2+x2^2)+(t3-t2)^4*v^4

A*x^2+B*x+C*y*x+D*y^2+E*y+F=0 　・・・式19
G*x^2+H*x+I*y*x+J*y^2+K*y+L=0　 　・・・式20
M*x^2+N*x+O*y*x+P*y^2+R*y+S=0 　・・・式21
　D、J、Pがいずれも０でないとき
　 y^2 の項を消去する。
J*(A*x^2+B*x+C*y*x+D*y^2+E*y+F)-D*(G*x^2+H*x+I*y*x+J*y^2+K*y+L)
=
(+A*J-D*G)*x^2+(B*J-D*H)*x+(C*J-D*I)*x*y+(E*J-D*K)*y+(F*J-D*L) ＝0 ・・・式22
　同様に
P*(A*x^2+B*x+C*y*x+D*y^2+E*y+F)-D*(M*x^2+N*x+O*y*x+P*y^2+R*y+S)
=
+(A*P-D*M)*x^2+(B*P-D*N)*x+(C*P-D*O)*x*y+(-D*R+E*P)*y+(-D*S+F*P)=0 ・・・式23

　式22、式23をそれぞれyの式にすると
(+A*J-D*G)*x^2+(B*J-D*H)*x+(C*J-D*I)*x*y+(E*J-D*K)*y+(F*J-D*L)=0
→y=-((A*J-D*G)*x^2+(B*J-D*H)*x-D*L+F*J)/((C*J-D*I)*x-D*K+E*J)　・・・式24

+(A*P-D*M)*x^2+(B*P-D*N)*x+(C*P-D*O)*x*y+(-D*R+E*P)*y+(-D*S+F*P)=0
→y=-((A*P-D*M)*x^2+(B*P-D*N)*x-D*S+F*P)/((C*P-D*O)*x-D*R+E*P)　・・・式25

式24−式25において分母を払うために
((C*P-D*O)*x-D*R+E*P)* ((C*J-D*I)*x-D*K+E*J)* (式24-式25)=0
とする

((C*P-D*O)*x-D*R+E*P)*(-((A*J-D*G)*x^2+(B*J-D*H)*x-D*L+F*J))
-((C*J-D*I)*x-D*K+E*J)*(-((A*P-D*M)*x^2+(B*P-D*N)*x-D*S+F*P))=0

これは
((C*D*G-A*D*I)*P+(A*D*J-D^2*G)*O+(D^2*I-C*D*J)*M) *x^3
+((A*D*J-D^2*G)*R+(-A*D*K-B*D*I+C*D*H+D*E*G)*P+(B*D*J-D^2*H)*O+(D^2*I-C*D*J)*N+(D^2*K-D*E*J)*M) *x^2
+((D^2*I-C*D*J)*S+(B*D*J-D^2*H)*R+(C*D*L-B*D*K-D*F*I+D*E*H)*P+(D*F*J-D^2*L)*O+(D^2*K-D*E*J)*N) *x
+(D^2*K-D*E*J)*S+(D*F*J-D^2*L)*R+(D*E*L-D*F*K)*P
=0

ここで
a= ((C*D*G-A*D*I)*P+(A*D*J-D^2*G)*O+(D^2*I-C*D*J)*M)
b= +((A*D*J-D^2*G)*R+(-A*D*K-B*D*I+C*D*H+D*E*G)*P+(B*D*J-D^2*H)*O+(D^2*I-C*D*J)*N+(D^2*K-D*E*J)*M)
c= +((D^2*I-C*D*J)*S+(B*D*J-D^2*H)*R+(C*D*L-B*D*K-D*F*I+D*E*H)*P+(D*F*J-D^2*L)*O+(D^2*K-D*E*J)*N)
d= +(D^2*K-D*E*J)*S+(D*F*J-D^2*L)*R+(D*E*L-D*F*K)*P
と置くと与式は、

　a*x^3 ＋b*x^2 ＋c*x ＋d＝０

と書ける。

→式の選び方で、「２次方程式の解を求める」問題になるのか、「３次方程式の解を求めざるを得ない」問題になるのか、ということになります。

２次方程式の一般解の公式は高校で習いました。

　a*x^2 ＋b*x ＋c ＝０

の一般解の公式は

x = （-ｂ+SQRT（b^2-4*a*c))/(2*a)
　or
x = （-ｂ−SQRT（b^2-4*a*c))/(2*a)

なので、３次方程式の一般解もたいしたことないだろうと想像しますが・・・
ちなみに、

　a*x^3 ＋b*x^2 ＋c*x ＋d＝０

の一般解は、

　x＝−【c*((−32^(4/3))*a*β^(2/3)−(92^(2/3))*a*b*β^(1/3)−18*a*b^2)+(2^(4/3))*(b^2)*(β^(2/3)+d*((272^(2/3))*(a^2)*β^(1/3)+54*(a^2)*b)−(2^(2/3))*(3^(3/2))*a*α*β^(1/3)+(2^(5/3))*(b^3)*(β^(1/3))−(23^(3/2))*a*b*α+4*b^4】/【−(23^(5/2))*(a^2)*α+163*(a^3)*d−54*(a^2)*b*c+12*a*b^3】　・・・式26

　x＝−【c*(a*((32^(4/3))*β^(2/3)−(2^(4/3))*(3^(3/2))*%i*β^(2/3))+a*b*(+(2^(2/3))*(3^(5/2))*%i*β^(1/3)+(92^(2/3))*β^(1/3))−36*a*b^2)+(b^2)*（+(2^(4/3))*(3^(1/2))*β^(2/3)−(2^(4/3))*β^(2/3))+d*((a^2)*(−(2^(2/3))*(3^(7/2))*%i*β^(1/3)−(272^(2/3)*β^(1/3))+108*(a^2)*b)+a*(+(92^(2/3)*%i*α*β^(1/3)+(2^(2/3))*(3^(3/2))*α*β^(1/3))+(b^3)*(−(2^(5/3))*(3^(1/2))* %i*β^(1/3)−(2^(5/3))*β^(1/3))−(43^(3/2))*a*b*α+8*b^4】/【−(43^(5/2))*(a^2)*α+324*(a^3)*d−108*(a^2)*b*c+24*a*b^3】　・・・式27

　x＝−【c*(a*((32^(4/3))*β^(2/3)+(2^(4/3))*(3^(3/2))*%i*β^(2/3))+a*b*(−(2^(2/3))*(3^(5/2))*%i*β^(1/3)+(92^(2/3))*β^(1/3)) −36*a*b^2)+(b^2)*（−(2^(4/3))*(3^(1/2))*β^(2/3)−(2^(4/3))*β^(2/3))+d*((a^2)*(+(2^(2/3))*(3^(7/2))*%i*β^(1/3)−(272^(2/3))*β^(1/3))+108*(a^2)*b)+a*(−(92^(2/3))*%i*α*β^(1/3)+(2^(2/3))*(3^(3/2))*α*β^(1/3))+(b^3)*(+(2^(5/3))*(3^(1/2))*%i*β^(1/3)−(2^(5/3))*β^(1/3))−(43^(3/2))*a*b*α+8*b^4】/【−(43^(5/2))*(a^2)*α+324*(a^3)*d−108*(a^2)*b*c+24*a*b^3】　・・・式28

見やすい形式は

数値計算をしたいとは思わないです。

数式を展開したり因数分解をしたりする操作を手作業で行っているとどうしても間違いが生じやすいのです。
そこで、数式処理ソフトを活用すると、この操作の煩わしさから解放されて、その上間違いが少なくなります。
数式処理ソフトで操作後に自分好みに加工する際に間違いが生じるという、情けないはなしなのですが・・。

無料の数式処理ソフトがインターネットから入手できます。

→「数式処理ソフトMaxima」

使い方の入門書としては、ブルーバックス（竹内薫著）があります。

その先をやってみたいと思ったら、「はじめてのMaxima（横田博史著）」があります。

お試しあれ。

音源定位の式の展開及び因数分解にMaximaをおおいに活用しました。
また3次方程式の一般解はMaximaを使って求めました。
Maximaでは4次方程式の一般解は求められません。（反応しない？）

一度、Mathematica（有料の数式処理ソフト）を使える機会があったときに、4次方程式の一般解を求めてみました。
Ａ４の紙に数ページにわたって印刷しなければ書ききれなきくらい長い解の数式が出てきていました、
「この根号はどこまで被っているのか？」というくらい長いし、また根号の中にさらに根号があり、読み解くのが大変です。

5次以上の方程式の一般解は存在しないそうです。
正確には
「代数学の基本定理によれば、任意の複素数係数方程式は複素数の中に根が存在するが、しかしながら5次以上の方程式には一般には代数的解法は必ずしも存在しない。
すなわち、一般の五次方程式に対して代数的な根の公式は存在しない。
もう少し詳しく書くと、5次の一般方程式の根を、その式の各項の係数と有理数の、有限回の四則演算及び有限回の根号をとる操作の組み合わせで表示することはできない。」
（Wikipedia）

ここらあたりは「ガロア理論」の領域ですね。

ページのトップへ

----- ここまで　20180429 -----

HOMEに戻る