3次元上の座標計算がうまく解けなくて困っています。

質問

3次元上の座標計算がうまく解けなくて困っています。

3次元上に2つの座標P(1)とP(2)があり、それぞれ単位ベクトルV(1)とV(2)を持っています。
P(1)とP(2)はそれぞれそのベクトル上に座標P(a)P(b)を持ちます。
つまり
P(a)=P(1)+sV(1)
P(b)=P(2)+tV(2)
s>0, t>0
と言う関係になります。
P(1)P(a)P(b)P(2)の順で線を結んだとき、P(1)P(a)P(b)がなす面上に半径Rの円を描きます。この円はP(1)P(a)が接線となりP(a)P(b)も接線となります。またP(1)P(a)での接点はP(1)になります。
同様にP(a)P(b)P(b)がなす面上に同じ半径Rの円を描きます。この円もP(b)P(2)とP(a)P(b)が接線となり、P(b)P(2)での接点はP(2)になります(図)。
上記の様になるsとtを求めたいのですが、どうすれば良いでしょうか？

一応、考えたのはP(1)P(a)P(b)でなす角度をθ、P(a)P(b)P(2)でなす角度をφとしたときに
cos(θ/2)=s/√(s^2+R^2)
sin(θ/2)=R/√(s^2+R^2)
から加法定理より
cos(θ)=(s^2-R^2)/(s^2+R^2)
同様に
cos(φ)=(t^2-R^2)/(t^2+R^2)

またP(b)P(a)への単位ベクトルをV(a)としたとき、
V(a)=(P(a)-P(b))/|P(a)-P(b)|
cos(θ)=V(1)・V(a)
同様に
cos(φ)=V(2)・-V(a)
となりますので、
(s^2-R^2)/(s^2+R^2)=V(1)・V(a)
(t^2-R^2)/(t^2+R^2)=V(2)・-V(a)
から解けないかと思ったのですが、上手く展開できませんでした。

よろしくお願いします。

回答No.4

stomachman · ベストアンサー

ANo.3 が字化けしちゃった…

"•"ってのは内積の記号"・"です。化けてるのはこの記号だけなんで、エディタで置き換えてナニしてください。すいません。

お礼

回答、ありがとうございます。
数学にうといので解析的、数値的の意味を勘違いして、返答してしまいました。
解析的に解くのが無理なら何か収束させていくような方法を取るしかないかなと思っていたので、参考になりました。

回答No.3

stomachman

ANo.2のコメントについてです。
数値的に答が出ればいいのなら、諦めることはありません。

まずは
[1] Pb(=P2+tV2)が与えられたとき、sは幾らになるか。
[2] Pa(=P1+sV1)が与えられたとき、tは幾らになるか。
の2本の方程式を連立してみました。その結果、解析的に簡単に解けるとは思えないんですが、数値計算なら難しくないようです。

まず、[1]を導出してみます。Pbが与えられると、Pb, P1,V1で平面が決まります。この平面に直交座標系を作りましょう。y軸の単位ベクトルとしてV1を使います。x軸の単位ベクトルをU1と書くことにすると、u1はV1と直交しているから
U1 = (b1-(b1•V1)V1)/|b1-(b1•V1)V1|
ただし•は内積で、また、
b1 = Pb-P1
です。
　円の中心C1を原点にしましょう。従って、点P1の座標は(±R,0)です。ご質問の図で見ると、どうも(-R,0)だけ考えればいいみたいですが、一応、P1=(-n1R,0) （n1∈{-1,1})、ということにします。（ベクトルとして書けばP1=C1-nRuです。）
点Pbの座標は(f,g)、すなわち
b1 = f1U1+g1V1
g1=b1•V1
f1=b1•U1-n1R
U1 = (b1-g1V1)/|b1-g1V1|
と表せます。b1•V1≧0, b1•U1≧0ですが、題意からするとb1•U1＞0ですね。

　さて、(f1,g1)を通る傾きαの直線L(α)
y = α(x-f1)+g1
と原点との距離√(x^2+y^2) はxが
x[0]=(αf1-g1)α/(1+α^2)
のとき最小になり、その最小値がRに等しいのだから、
R^2 = x[0]^2+(α(x[0]-f1)+g)^2
これはαに関する2次方程式で、解は
α = (g1±R)/f1
です。円には接線が2本ひけるから、コタエも2個出てくる。複号だとややこしいんで、
k1∈{-1,1}
ということにして
α = (g1+k1R)/f1
と書きましょう。この直線L(α)が直線
x=-n1R
と交わる点がPaであり、そのy座標がs、つまり
s = α(-n1R-f1)+g1
なので、
s= (k1- n1(g1+k1R))R/f1

　ここに出てくるf1, g1をP1, P2, t , V1, V2, Rで表せばいい。
b1 = P2-P1+tV2
だから、
g1 = b1•V1 = (P2-P1)•V1 + t (V1•V2)
f1=b1•U1-n1R　= b1•(b1-g1V1)/|b1-g1V1|-n1R
ここに
|b1-g1V1|=√((b1-g1V1)•(b1-g1V1))
です。これでsをtで表した式ができた訳です。

　この式に出てくる(P1, P2,V1, V2, s, t, n1, k1)をそれぞれ(P2, P1, V2, V1, t, s, n2, k2)に置き換えれば、tをsで表す式ができ、この２本で連立方程式になっています。（が、平方根のせいで、解析的に解くのは難しいっぽい。）

　なら数値的に解く。やり方はいろいろありますけど、まずは一番楽ちんな手として、適当な出発値、たとえばt=1から始めて、sを計算し、それを使ってtを計算し、... とやってみる。もちろん、これで発散しちゃうようであればもっと凝ったやり方を使う必要がありますが。
　n1, n2, k1, k2の選び方で2^4通りの組み合わせがありますんで、まずはそれぞれやってみるですかね。

　という訳で、excelで試しました。P1, P2, V1, V2, n1, n2, k1, k2, Rを定数として表で与えておいて、１行で、
t → ベクトルb1 → b1•V1→ベクトルb1-(b1•V1)→|b1-(b1•V1)|→ベクトルU1→g1, f1→ s → ベクトルb2 → b2•V2→ベクトルb2-(b2•V2)→|b2-(b2•V2)|→ベクトルU2→g2, f2→ t
を計算し、これを次の行の先頭に持ってくる。以下同様。
　定数はデタラメに
P1=(0,2,0),P2=(5,3,7), V1=(1,0,0), V2=(0,(√2)/2,(√2)/2), R=1.3
と与えて、試しに
n1=1, k1=1, n2=-1, k2=-1 としたところ、t, sは
1.0000, -0.8231
0.4573, -0.7417
0.3998, -0.7242
0.3861, -0.7195
0.3824, -0.7182
0.3813, -0.7178
0.3810, -0.7177
0.3809, -0.7177
0.3809, -0.7176
0.3809, -0.7176
:
と収束してくれました。イケるかもです。ただしstomachmanは計算間違いの常習犯ですんで、しっかりチェックして下さいな。

回答No.2

stomachman

滑車を使ってワイヤーを張る、みたいな話かな？
　試しに三角関数なしでやってみましたら、s,tの連立方程式は作れたけど、どうも解ける気がしません。数値計算ならなんとかなるでしょう。

お礼

そうですね。ワイヤを張るみたいな感じです。
そういう数値計算のプログラムを作りたかったんですが、やはり難しそうですね。

回答No.1

nag0720

P(1),P(a),P(b),P(2)は同一平面上にあるとしていいのですか？

P(1)からのベクトルV(1)方向と、P(2)からのベクトルV(2)方向がねじれていると同一平面にはならないですが・・・

補足

同一平面でない(V(1)とV(2)がねじれている)場合を求めたいのです。
同一平面の場合は何とか自力で出来たので、ねじれのケースでお願いします。