5元連立方程式を解きたい

質問

複雑な方程式を解きたいと考えています。

方程式は、未知数は5つ、関係式も5つあるので、
解析的に解くことが可能な5元連立方程式だと思います。
複雑、といっても、高度な関数などが登場するのではなく
(せいぜい三角関数程度)、項の数がかなり多いだけです。

とりあえず手計算を試みたのですが、
項を書き下すことすら難しいような状況で、
そのまま計算を進めても、結局、解まで辿り着けたことはありません
(数人で挑戦したのですが、みんな駄目でした)。

フリーの数式処理ソフト Maxima も導入して計算させてみたのですが、
エラーを吐いてしまい、やはり解は求まりませんでした
(誤植等が無いか何度も確認しましたが、駄目でした)。

他の、有償のソフトならば解ける可能性はあるでしょうか?
また、方程式が解析的に解けるものなのかどうかを
判断する方法などがありましたら、お教え下さい。

回答No.5

inara · ベストアンサー

pdfファイルは、URLを直接アドレスバーに入力することで読めました。最小２乗法ですね。私も仕事で、測定データを非線形関数でfittingするときに同じような計算（残差２乗和を偏微分）をしたことがあります。

数式処理ソフト（Maple）で挑戦してみました。
dSax～dSby までは計算可能です（ものすごい長い式）が、この連立方程式の解は出てきませんでした。
そこで、手動操作で変数を減らすために、dSax=0 を a2x について解かせてみましたが、この段階ですでに解が出てきません。

ここであきらめるのも癪（しゃく）なので、途中の変数の次数を調べてみました。Mapleには式の中の変数の最高次数を求めるコマンド degree(式, 変数）があります。

ある変数（文字式）を有理文字式（文字式/文字式）で表したとき、そこに含まれている５変数 a2x, a2y, a2z, b2x, b2y の最高次数をそれぞれ N1, N2, N3, N4, N5, N6 として、その変数を 分子（N1,N2,N3,N4,N5) /分母（N1,N2,N3,N4,N5) と書き表すことにします。
すると、L1～L6 は 分子(4,4,4,2,2)/分母(4,4,4,0,0)
でしたので、その次の計算
S = (L1-d1)^2 + (L2-d2)^2 + (L3-d3)^2 + (L4-d4)^2 + (L5-d5)^2 + (L6-d6)^2
では、S は分子(8,8,8,4,4)/分母(8,8,8,0,0) となるはずです（Mapleでは、この段階で最高次数計算が不能となった）。

さらに、S を微分したものは、商の微分法を考えれば、分子、分母とも、最高次数は 16 になると思われます（Mapleでは文字式自身は出ますが、式が複雑すぎて、項別にまとめたり、因数分解はできません。展開された恐ろしく長い式を手計算でまとめるのはちょっと無理）。

したがって問題の５元連立方程式は、5個の「16次多項式/16次多項式=0」を解くということになります。特殊な場合には解析解があるかもしれませんが、一般的には解析的に解けないかもしれません。

Mapleは文字式操作は得意かと思っていましたが、これくらいで計算不能となるとは思いませんでした。私の使ったMapleがVersion5で、しかもStudent版（学生用に機能制限された廉価版）だったせいかもしれませんが、もし、Mapleの最新の正式版をお持ちの方があれば、以下にMapleのコードを書きますので、解は出るのか、どこまで文字式計算可能か調べてみてください。

alpha1:=a1x1*a2x + a1y1*a2y + a1z1*a2z:beta1:=a1x1*(b2x-b1x1) + a1y1*(b2y-b1y1) + a1z1*(b2z-b1z1):gamma1:=a2x* (b2x-b1x1) + a2y* (b2y-b1y1) + a2z* (b2z-b1z1):alpha2:=a1x2*a2x + a1y2*a2y + a1z2*a2z:beta2:=a1x2*(b2x-b1x2) + a1y2*(b2y-b1y2) + a1z2*(b2z-b1z2):gamma2:=a2x* (b2x-b1x2) + a2y* (b2y-b1y2) + a2z* (b2z-b1z2):alpha3:=a1x3*a2x + a1y3*a2y + a1z3*a2z:beta3:=a1x3*(b2x-b1x3) + a1y3*(b2y-b1y3) + a1z3*(b2z-b1z3):gamma3:=a2x* (b2x-b1x3) + a2y* (b2y-b1y3) + a2z* (b2z-b1z3):alpha4:=a1x4*a2x + a1y4*a2y + a1z4*a2z:beta4:=a1x4*(b2x-b1x4) + a1y4*(b2y-b1y4) + a1z4*(b2z-b1z4):gamma4:=a2x* (b2x-b1x4) + a2y* (b2y-b1y4) + a2z* (b2z-b1z4):alpha5:=a1x5*a2x + a1y5*a2y + a1z5*a2z:beta5:=a1x5*(b2x-b1x5) + a1y5*(b2y-b1y5) + a1z5*(b2z-b1z5):gamma5:=a2x* (b2x-b1x5) + a2y* (b2y-b1y5) + a2z* (b2z-b1z5):alpha6:=a1x6*a2x + a1y6*a2y + a1z6*a2z:beta6:=a1x6*(b2x-b1x6) + a1y6*(b2y-b1y6) + a1z6*(b2z-b1z6):gamma6:=a2x* (b2x-b1x6) + a2y* (b2y-b1y6) + a2z* (b2z-b1z6):

t11:= (beta1-alpha1*gamma1)/(1-alpha1^2):t21:= (alpha1*beta1-gamma1)/(1-alpha1^2):t12:= (beta2-alpha2*gamma2)/(2-alpha2^2):t22:= (alpha2*beta2-gamma2)/(2-alpha2^2):t13:= (beta3-alpha3*gamma3)/(3-alpha3^2):t23:= (alpha3*beta3-gamma3)/(3-alpha3^2):t14:= (beta4-alpha4*gamma4)/(4-alpha4^2):t24:= (alpha4*beta4-gamma4)/(4-alpha4^2):t15:= (beta5-alpha5*gamma5)/(5-alpha5^2):t25:= (alpha5*beta5-gamma5)/(5-alpha5^2):t16:= (beta6-alpha6*gamma6)/(6-alpha6^2):t26:= (alpha6*beta6-gamma6)/(6-alpha6^2):r1x1:= t11*a1x1 + b1x1:r1y1:= t11*a1y1 + b1y1:r1z1:= t11*a1z1 + b1z1:r2x1:= t21*a2x + b2x:r2y1:= t21*a2y + b2y:r2z1:= t21*a2z + b2z:r1x2:= t12*a1x2 + b1x2:r1y2:= t12*a1y2 + b1y2:r1z2:= t12*a1z2 + b1z2:r2x2:= t22*a2x + b2x:r2y2:= t22*a2y + b2y:r2z2:= t22*a2z + b2z:r1x3:= t13*a1x3 + b1x3:r1y3:= t13*a1y3 + b1y3:r1z3:= t13*a1z3 + b1z3:r2x3:= t23*a2x + b2x:r2y3:= t23*a2y + b2y:r2z3:= t23*a2z + b2z:r1x4:= t14*a1x4 + b1x4:r1y4:= t14*a1y4 + b1y4:r1z4:= t14*a1z4 + b1z4:r2x4:= t24*a2x + b2x:r2y4:= t24*a2y + b2y:r2z4:= t24*a2z + b2z:r1x5:= t15*a1x5 + b1x5:r1y5:= t15*a1y5 + b1y5:r1z5:= t15*a1z5 + b1z5:r2x5:= t25*a2x + b2x:r2y5:= t25*a2y + b2y:r2z5:= t25*a2z + b2z:r1x6:= t16*a1x6 + b1x6:r1y6:= t16*a1y6 + b1y6:r1z6:= t16*a1z6 + b1z6:r2x6:= t26*a2x + b2x:r2y6:= t26*a2y + b2y:r2z6:= t26*a2z + b2z:

L1:= (r1x1-r2x1)^2 + (r1y1-r2y1)^2 + (r1z1-r2z1)^2:L2:= (r1x2-r2x2)^2 + (r1y2-r2y2)^2 + (r1z2-r2z2)^2:L3:= (r1x3-r2x3)^2 + (r1y3-r2y3)^2 + (r1z3-r2z3)^2:L4:= (r1x4-r2x4)^2 + (r1y4-r2y4)^2 + (r1z4-r2z4)^2:L5:= (r1x5-r2x5)^2 + (r1y5-r2y5)^2 + (r1z5-r2z5)^2:L6:= (r1x6-r2x6)^2 + (r1y6-r2y6)^2 + (r1z6-r2z6)^2:

S:=(L1-d1)^2+(L2-d2)^2+(L3-d3)^2+(L4-d4)^2+(L5-d5)^2+(L6-d6)^2:

dSax:=diff(S, a2x):dSay:=diff(S, a2y):dSaz:=diff(S, a2z):dSbx:=diff(S, b2x):dSby:=diff(S, b2y):

solve({dSax=0, dSay=0, dSaz=0, dSbx=0, dSby=0}, {a2x, a2y, a2z, b2x, b2y});

お礼

自分が Maxima で計算してみたときも、
dSax~dSby, つまり微分までは出来たのですが、
これを展開させようとすると詰まってしまう状況でした
(もちろん、方程式を計算させても解は出ない)。

しかしまぁ、なんと最高次数 16 ですか!
5 程度かなとか考えていた自分は甘過ぎたようです。
解析解を望むのは、ちょっと難しそうですね。

正式版の Maple を入手できたら、書いていただいたコードを
計算させてやりたいと思います。

計算の挑戦のみならず、丁寧な考察までしていただいて、
本当にありがとうございました。

回答No.4

inara

equation.pdf のほうは、「このページは準備中です」と出て読めないのですが、これは equation.txt と同じでしょうか。 equation.txt は読めます（しかし大変な5元連立方程式ですね）。

関数の定義の部分は、Mapleと定義式が若干違いますが理解可能です。diff(S(), a2x, 1) の意味は、関数 S() を a2x で1階微分するという意味ですね？　solve のところは括弧の使い方の違いだけですので、これも理解できます。

コードの意味は十分理解できますので、さっそく挑戦してみます。質問はこのまま開いておいてください。

補足

うーん、URLは確実に合ってるんですが、なぜかクリックするのではつながらないようです。
このURLを直接アドレスバーに入力していただければ、OKみたいです。
http://www.geocities.jp/zacolih0083/study/equation.pdf
pdfファイルでは、そもそもこの方程式は一体何なのか、ということを解説しています
(仲間内で公開したものなので、解説不足な点があるとは思いますが
　数式がtxtより読みやすくなってます)。

関数定義は、その通りです。
pdfが読めれば、確認できると思います。

ありがとうございます。よろしくお願いします。

回答No.3

inara

計算過程でなく、結果だけが必要なのであれば、手元に数式処理ソフト（Maple）があるので、5元連立方程式を教えていただければ、数式解を解いてみましょう。方程式によっては解析解は見つからないかもしれませんが。

補足

では、お言葉に甘えて、方程式を紹介させてもらいます。
http://www.geocities.jp/zacolih0083/study/equation.pdf
「これを解けば良い」という形にまとまっていないので、
分かりにくいと思います。すみません。
Maximaで計算させたコードも載せておきます。
http://www.geocities.jp/zacolih0083/study/equation.txt
おそらくMapleでは文法が違うので、そのままでは使えないとは思いますが、
参考になればなりよりです。
お時間があるときで構いませんから、
チャレンジしてみていただけると大変助かります。

回答No.2

麻野 なぎ（@AsanoNagi）

近似的な数値解でよれば、（かつ、Excel がお手元にあれば）Excel のソルバーをつかうというのもひとつの手かと思います。
下記の URL に少し触れてあるのを見つけました。
（ページの下の方の、「大人は誰も教えてくれない連立方程式の解法」）

ただ、数値解析なので、初期条件によってはうまく収束しないとか、そういったことも起こることがあります。

参考URL：

http://www.sun-inet.or.jp/~yaneurao/rsp/rsp29to2F.html

補足

実は、数値計算はすでに行なっているのです。
一応、C言語を扱えるので、それなりに工夫して近似解を
求められるようにはなったのですが、
究極的にやりたいことは、方程式の各種パラメータを様々に変更し、
その度に解を求める、というものなのです。
その数も数千、数万という数なので、
数値計算ではなかなか実用的なほど効率化できないので、
解析的に求められないか、ご相談させてもらいました。

回答No.1

uninin

次数はいかほどでしょうか。

補足

きちんと把握できていないのですが、
おそらく 5 程度だと思います。