3次方程式のカルダノ（Cardano)法によるエクセルでの計算方法について

皆さんよろしくお願いいたします。
標題の計算をエクセルで行ってみました。
ところが、計算結果が係数によって合うこともあるのですが、
3つの異なる複素解が出てしまうことが有り、ありえない数値が出ることで困っています。
具体的には、3次方程式をax^3+bx^2+cx+d=0としたとき係数に次の数値を入れたときに起こります。
a=1/120、b=1/10,c=1/2,d=1
同じようなご経験のある方、ご教示いただければ幸いです。
計算方法は記号など以下のサイトを引用しました。
http://www2.biglobe.ne.jp/~ytajima/cardanos_form …
計算内容は以下のとおりです。
A1,B1,C1,D1セルにそれぞれの係数a,b,c,dを代入します。
A2=B1/A1---p=b/aの計算
A3=C1/A1---q=c/aの計算
A4=D1/A1---r=d/aの計算
A5=-(A2^2)/3+A3---m=-p^2/3+qの計算
A6=(2*(A2^3)/27)-(A2*A3/3)+A4---n=2p^3/27-pq/3+rの計算
A7=IMPOWER(IMSUM((-A6/2),IMPOWER((A6^2/4+A5^3/27),1/2)),1/3)---u=(-n/2+(n^2/4+m^3/27)^0.5)^(1/3)の計算
A8=IMPOWER(IMSUB((-A6/2),IMPOWER((A6^2/4+A5^3/27),1/2)),1/3)---v=(-n/2-(n^2/4+m^3/27)^0.5)^(1/3)の計算
A9=COMPLEX(A10,A11)---ω=(-1+√3i)/2の計算
A10=-1/2---Re(ω)
A11=SQRT(3)/2---Im(ω)
A12=IMPRODUCT(A9,A9)---ω^2の計算
A13=IMSUB(IMSUM(A7,A8),(A2/3))---解1=u+v-p/3の結果
A14=IMSUB(IMSUM(IMPRODUCT(A9,A7),IMPRODUCT(A8,A12)),(A2/3))---解2=ωu+vω^2-p/3の結果
A15=IMSUB(IMSUM(IMPRODUCT(A9,A8),IMPRODUCT(A7,A12)),(A2/3))---解3=uω^2+vω-p/3の結果

No.7 ベストアンサー 回答者： Tacosan 回答日時： 2008/04/14 23:38

あれ? 手元の Excel 2007 SP1 (#2 と同じ環境です) で試したんですが, 期待通りの結果が出てますよ.

例えば
A1=-1/120
B1=1/10
C1=-1/2
D1=1
A2=B1/A1
A3=C1/A1
A4=D1/A1
A5=-(A2^2)/3+A3
A6=(2*(A2^3)/27)-(A2*A3/3)+A4
A7=IMPOWER(IMSUM((-A6/2),IMPOWER((A6^2/4+A5^3/27),1/2)),1/3)+A7:A69
A8=IMPRODUCT(-1/3,IMDIV(A5,A7))
A9=COMPLEX(A10,A11)
A10=-1/2
A11=SQRT(3)/2
A12=IMPRODUCT(A9,A9)
A13=IMSUB(IMSUM(A7,A8),(A2/3))
A14=IMSUB(IMSUM(IMPRODUCT(A9,A7),IMPRODUCT(A8,A12)),(A2/3))
A15=IMSUB(IMSUM(IMPRODUCT(A9,A8),IMPRODUCT(A7,A12)),(A2/3))
とすると (つまり 4 ですね)
A13=4.64437070925218
A14=3.67781464537391+3.50876191956744i
A15=3.67781464537391-3.50876191956744i
が得られ, それぞれの残差は
2.88657986402541E-15
-6.21724893790088E-15-3.10862446895044E-15i
-6.21724893790088E-15+3.10862446895044E-15i
です.

この回答へのお礼

ご回答頂きありがとうございます。
どうやら、小生のエクセルに計算違いや入力ミスがあったようです。
おかげさまで間違いに気づくことが出来ました。
ありがとうございます。

お礼日時：2008/04/15 00:34

No.6 回答者： noname#101087 回答日時： 2008/04/14 20:23

>u=(-n/2+(n^2/4+m^3/27)^0.5)^(1/3)

>v=(-n/2-(n^2/4+m^3/27)^0.5)^(1/3)

この {u, v} のペアは
　　・いずれも実数
　　・互いに共役な(非実)複素数
のいずれか。

どちらにせよ、このペアが実零点(x = u+v-p/3)の一つを与えるわけです。
　

No.5 回答者： Tacosan 回答日時： 2008/04/14 13:08

今は手元に計算できる環境がないので, 回答ではなく指摘だけ:

IMPOWER で 3乗根を計算するというのは, 数学的には
z = r e^(iθ) から z^(1/3) = r^(1/3) e^(iθ/3) を計算する
ことと等しいはずです... というか結果から見る限りこの計算をしているようです.
で, z として負の実数, 例えば -1 を入れると
-1 = 1 ・ e^(iπ) だから z^(1/3) = (-1)^(1/3) として
1^(1/3) ・ e^(iπ/3) = 1/2 + (√3/2)i
が得られます. これは最初の例でも本質的に同じで,
1.17398499670533+2.03340166161721i
の実部と虚部の比は 1:√3 です.
これに対し正の実数を入れると 3乗根として (正の) 実数となり, ここで uv = -m/3 という関係を崩すことになります. だから v = -m/(3u) と計算する, と書きました.
で今はどのように計算しているのでしょうか. 全体を改めて書いていただいて, その上で例えば (2) に対して期待される結果と実際に得られた値を書いていただけるとありがたいです.

この回答へのお礼

ご回答頂きありがとうございます。
数学的な計算方法など懇切丁寧にご指導頂きありがとうございます。
今どのように計算しているかというと
v=(-n/2-(n^2/4+m^3/27)^0.5)^(1/3)の計算
A8=IMPOWER(IMSUB((-A6/2),IMPOWER((A6^2/4+A5^3/27),1/2)),1/3)
を
A8=IMPRODUCT(-1/3,IMDIV(A5,A7))
としてみました。これはuv = -m/3を満たすようにvを求めるようにしたものです。しかしながら結果は求めたい値に対し以下のようになりました。
1)a=1/120,b=1/10,c=1/2,d=1（解決済み）
　自分の計算結果と得たい結果は精度の問題はありますが一致しました。
　x1=-4.6443707092516、X2=-3.67781464537391+3.50876191956745i
　x3=-3.67781464537391-3.50876191956745i
2)a=8/120,b=4/10,c=1,d=1
　自分の計算x1=0,x2=-3+2.44948974278319i,x3=-3-2.44948974278319i
　得たい結果x1=-2.3221853357179962,
　　　　　　x2=-1.8389073171410022+1.754380971381805i
　　　　　　x3=-1.8389073171410022-1.754380971381805i
3)a=125/120,b=25/10,c=5/2,d=1
　自分の計算x1=5.107025913275752E-015,
　　　　　　x2=-1.2+0.979795897113268i
　　　　　　x3=-1.2-0.979795897113268i
　得たい結果x1=-0.9288741413663972,
　　　　　　x2=-0.7355629289328063+0.7017523842103995i,
　　　　　　x3=-0.7355629289328063-0.7017523842103995i
4)a=-1/120,b=1/10,c=-1/2,d=1
　自分の計算x1=0,x2=6+4.89897948556636i,x3=6-4.89897948556636i
　得たい結果x1=4.644371011781658
　　　　　　x2=3.677814734101803+3.508761733998084i
　　　　　　x3=3.677814734101803-3.508761733998084i
5)a=-8/120,b=4/10,c=-1,d=1
　自分の計算x1=0,x2=3+2.44948974278319i,x3=3-2.44948974278319i
　得たい結果x1=2.3221853357179962,
　　　　　　x2=1.8389073171410022+1.754380971381805i
　　　　　　x3=1.8389073171410022-1.754380971381805i
6)a=-125/120,b=25/10,c=-5/2,d=1
　自分の計算x1=-5.10702591327872E-015
　　　　　　x2=1.2+0.979795897113268i
　　　　　　x3=1.2-0.979795897113268i
　得たい結果x1=0.9288741413663872
　　　　　　x2=0.7355629289328063+0.7017523842103995i
　　　　　　x3=0.7355629289328063-0.7017523842103995i

お礼日時：2008/04/14 14:16

No.4 回答者： noname#101087 回答日時： 2008/04/13 20:01

タイプミス。

　　v^3 = 4*{-1 + SQRT(5)}
　　w^3 = 4*{-1 - SQRT(5)}

No.3 回答者： noname#101087 回答日時： 2008/04/13 19:41

こちらの EXCEL は複素関数無しですけど.... 。

下記の問題なら実数演算で用を足せます。

-------------------------
原方程式：　x^3 + 12x^2 + 60x +120 = 0
x = y-4 として：　y^3 + 12y + 8 = 0
y = v + w として：
　　v^2 + w^2 + 8 =0
　　v*w = -12
この二式から：
　　v^6 + 8v^3 - 64 = 0
したがって：
　　v^3 = 4*{-1 + SQRT(5)}
　　v^3 = 4*{-1 - SQRT(5)}
あとは逆行。
　　y = -0.644...
　　x = -4.644...
…でした。

No.2 回答者： Tacosan 回答日時： 2008/04/13 18:35

あれ? 手元の Excel 2007 で計算した限りでは

A8=IMDIV(A5,IMPRODUCT(-3,A7))
としたら
A13 = -4.64437070925216
となり, 元の方程式に代入したときの残差は2.22044604925031E-15-2.21695285066186E-16i
になりました. これなら精度も十分でしょう.
ちなみにこの問題の場合 m = 12 で
u^3 = 4.94427190999916,
v^3 = -12.9442719099992
となっています. これらに対して IMPOWER で 3乗根を計算すると u の方は 1.70359928415849 と実数が得られますが, v は 1.17398499670533+2.03340166161721i
と複素数になり, uv = -m/3 を満たしません.

この回答への補足

お気づきの点『ちなみにこの問題の場合 m = 12 で～
u の方は 1.70359928415849 と実数が得られますが, v は 1.17398499670533+2.03340166161721i
と複素数になり, uv = -m/3 を満たしません.』
おっしゃるとおりですね、解は合っているのに制約条件を満たさないとは
不可思議に思います。自分も検索したり、エクセルの式を色々
いじくってますが、うまくいきません。

補足日時：2008/04/14 12:29

この回答へのお礼

ご回答頂きありがとうございます。
質問の係数a=1/120,b=1/10,c=1/2,d=1については、求めようとした解が出てきました。
ありがとうございます。
しかし、まだ以下の2)～5)の係数では、おかしな解が計算されてしまいます。
求めたい係数の種類は解決済みの1)の係数も含め、以下のとおりです。
1)a=1/120,b=1/10,c=1/2,d=1（解決済み）
2)a=8/120,b=4/10,c=1,d=1
3)a=125/120,b=25/10,c=5/2,d=1
4)a=-1/120,b=1/10,c=-1/2,d=1
5)a=-8/120,b=4/10,c=-1,d=1
6)a=-125/120,b=25/10,c=-5/2,d=1
結果を以下のHPで確認してみました。
http://homepage3.nifty.com/Architect/TE/ThreeEqu …
ちなみに上の1)～5)は制御工学の教科書に解が載っており、
上のサイトでの計算結果と一致します。
実用的に使用するために、どのような係数に対しても正しい結果を
エクセルで出力したいのですが、うまくいきません。
度々申し訳ありませんが、お気づきになられたことなど有りましたら
ご教示いただきたくお願いいたします。

お礼日時：2008/04/14 12:25

No.1 回答者： Tacosan 回答日時： 2008/04/13 00:08

Cardano の公式で 3次方程式を解く途中, u^3 と v^3 を求めてから u, v を計算するんですが, この u と v には

uv = -m/3
という制約が付いています (というかもともとこの式を満たさなければならない).
つまり, u と v を u^3, v^3 からそれぞれ独立に求めてはダメで, 一方を求めたあともう一方は uv = -m/3 から求めなきゃならないような気がします.

この回答へのお礼

ご回答頂きありがとうございます。
おっしゃるとおりuv = -m/3の制約があります。ご指摘ありがとうございます。
しかしなが、制約については、以下の計算手順を踏んで、解を計算しています。
u,vについては次式になります。
u=(-n/2+(n^2・4+m^3/27)^0.5)^(1/3)=Aとおく
v=(-n/2-(n^2・4+m^3/27)^0.5)^(1/3)=Bとおく
するとu,vの根はそれぞれ3つの解がでます。
ω=-1/2+(√3)i/2とすると
u=A,Aω,Aω^2
v=B,Bω,Bω^2
これらを組み合わせると9通りの解(u+v-p/3)ができますが、
ここでuv = -m/3を満たす解は次の3通りしか有りません。
u+v=A+B,Aω＋Bω^2,Aω^2+Bω
これは計算に反映されているので、ここは問題ないのではないかと
思われます。
ためしにuのみ求めてv=-m/(3u)として計算してみましたが、
解は合いませんでした。
せっかくアドバイスを頂いたのに申し訳有りませんが、
他にお気づきの点などありましたら、なんでも結構ですのでご教示
いただきたくお願いいたします。

お礼日時：2008/04/13 11:50