算数 数学

音源定位 update 20180429

音源定位って?

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

 

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

 

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

音源1個にマイク1個の場合

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

音源1個にマイク2個の場合

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

 

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

 

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

 

(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 -----

音源1個のマイク3個の場合

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

 

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

 

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

 

v : 音速(既知、気温により決まる)  

 

同様に、マイク2,3に音が到達する時間差(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次方程式からは通常解が2個得られるが、
   このx0をPに代入して 得られるt0について、
   t0 < t1, t2, t3 を満たさないものは除外する。

 

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

 

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

 

ページのトップへ

 

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

虚像を消す方法は?

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

 

そこで、マイクを6個準備して、3個づつ2個のグループでそれぞれ音源の位置を計算するとします。
今、マイクM1、M2、M3のグループで求められる音源のX座標を、x=b/a(真の音源のx座標)とx=d/c(虚像)、
マイクM4、M5、M6のグループで求められる音源のX座標を、x=b/a(真の音源のx座標)とx=-f/e(虚像)とします。

 

 

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

 

そこで、
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座標が導かれました。

 

Aは0でない実数のとき、 A (cx-d)= (ex-f)とする。いかなるxの値にたいしても、(Ac-e)x-(Ad-f)=0を
成り立たせるためには Ac=eかつAd=fであるから、A=e/c、A=f/d  ∴ e/c=f/d  → cf=de

 

 

マイクは6個もいるのか?

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

 

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

 

 

この場合、M1、M2、M3でひとつのグループ、M3、M4、M5でもうひとつのグループとして計算をします。
ひとつのグループ内のマイクの並び方は一直線です。

 

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

 

 

 

ページのトップへ

 

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

計算してみます

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

 

 

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

 

M1:(x,y)=(0m,1m)、測定時刻=12.700秒
M2:(x,y)=(3000m,-1m)、測定時刻=8.985秒
M3:(x,y)=(1000m,30m)、測定時刻=10.720秒
M4:(x,y)=(2000m,-30m)、測定時刻=9.5530秒

 

M1〜M2〜M3の組み合わせによる音源座標の計算は

 

 

及び

 

 

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

 

次にM1〜M2〜M4の組み合わせによる音源座標を計算します。

 

 

及び

 

 

このように計算結果を並べてみると、真の音源座標はなんとなく(x,y)=(3000m、3000m)と予想できます。

 

実際、音源座標は

 

 

に設定していました。

 

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

 

マイクの組み合わせM1〜M2〜M3による計算の地中で使用した係数 E1,F1,G1と、
M1〜M2〜M4での係数 E2,F2、G2を用いて

 

 

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

 

音源とマイクの座標から距離を計算して、それを音速で割った値を測定時刻に使うという、音源座標計算で誤算が出ないと思われる条件なのに、計算してみると、設定した音源座標と計算した音源座標を見比べてみると誤差が生じています。これは、各マイクでの時刻の測定精度を1000分の1からどんどん上げていくと、0になります。

 

ページのトップへ

 

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

 

気づきましたか?

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

 

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

 

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

 

です。
一見、音源が音波を発した時刻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
が導かれました。

 

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

 

まとめると、

ただし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がいずれも0でないとき
  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^3b*x^2c*xd=0

 

と書ける。

 

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

 

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

 

 a*x^2 +b*x +c =0

 

の一般解の公式は

 

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

 

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

 

 a*x^3 +b*x^2 +c*x +d=0

 

の一般解は、

 

 

 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次方程式の一般解を求めてみました。
A4の紙に数ページにわたって印刷しなければ書ききれなきくらい長い解の数式が出てきていました、
「この根号はどこまで被っているのか?」というくらい長いし、また根号の中にさらに根号があり、読み解くのが大変です。

 

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

 

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

 

ページのトップへ

 

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

 

HOMEに戻る