どこが一番近い？座標問題

質問

複数のポイントが座標上に在るとします。
例）
A　1,15
B　4,2
C　3,7
D　9,10
E　5,5

このような複数のポイントに一番合計距離が短くなる場所を求めるにはどう考えたらよいでしょうか？

回答No.11

178-tall

更なる蛇足。試したシナリオ。

距離和、
　g(x,y) = Σ gm(x,y)
　gm(x,y) = √{ (x-xm)^2 + (y-ym)^2 }　　: m = 1～5
にて、
　∂g(x,y)/∂x = Σ (x-xm)/gm(x,y)　…　gx
　∂^2g(x,y)/∂x^2 = Σ (y-ym)^2/{ gm(x,y) }^3　…　gxx
　∂^2g(x,y)/∂x∂y = ∂^2g(x,y)/∂y∂x
　= -Σ (x-xm)(y-ym)/{ gm(x,y) }^3　…　gxy = gyx
　∂g(x,y)/∂y = Σ (y-ym)/gm(x,y)　…　gy
　∂^2g(x,y)/∂y∂y = Σ (x-xm)^2/{ gm(x,y) }^3　…　gyy
を使っての逐次解法。

[ターゲット]　gx と gy が零なる (x, y) の探索。

[試点 (xt, yt) での次の試点 (xn, yn) ]　試点における {gx, gxx, gxy, gyy} を用いて、
　xn = xt - (gyy*gx - gxy*gy) / D
　yn = yt + (gxy*gx - gxx*gy) / D
　　: D = gxx*gyy - gxy^2

… (収束は遅め。ワントライで一桁の精度改良という感じ)

回答No.10

178-tall

蛇足をひとつ。

ANo.9 の「逐次解法」に使う gm の偏微分は、

　∂g(x,y)/∂x = Σ (x-xm)/gm(x,y)　…　gmx
　∂^2g(x,y)/∂x^2 = Σ (y-ym)^2/{ gm(x,y) }^3　…　gmxx
　∂^2g(x,y)/∂x∂y = ∂^2g(x,y)/∂y∂x
　= -Σ (x-xm)(y-ym)/{ gm(x,y) }^3　…　gmxy = gmyx
　∂g(x,y)/∂y = Σ (y-ym)/gm(x,y)　…　gmy
　∂^2g(x,y)/∂y∂y = Σ (x-xm)^2/{ gm(x,y) }^3　…　gmyy

らしい。

回答No.9

178-tall

ANo.6 の「フェルマ点証明の解析的方法」を眺めつつ、スプレッドシートにて「逐次解法」の単純繰り返しを試行。

まず、
　g(x,y) = Σ gm(x,y)
　gm(x,y) = √{ (x-xm)^2 + (y-ym)^2 }　　: m = 1～5
とする。

x, y にて偏微分すると、
　g'x(x,y) = Σ (x-xm)/gm(x,y)
　g''x(x,y) = Σ (y-ym)^2/{ gm(x,y) }^3
　g'y(x,y) = Σ (y-ym)/gm(x,y)
　g''y(x,y) = Σ (x-xm)^2/{ gm(x,y) }^3
らしい。

これらの勘定から、Newton 「逐次解法」で g'x(xo,yo) = g'y(xo,yo) = 0 に接近を試みると？

A　1,15、B　4,2、D　9,10　の「三角形」内のフェルマー点は、(6.072, 9.560) だった。
それを初期値とした Newton の「逐次解法」の 10 回目では (xo,yo) = (4.121, 6.568) 。

回答No.8

f272

#7です。

ちなみに#7に書いた答えはエクセルのソルバーで計算させた。
#4さんが試みたように2分法を使うのなら，こんな感じです。
Sub bisection(j, n, X, Y, minX, midX, maxX, minY0, midY0, maxY0, flag)
If flag = 1 Then
XX = minX
ElseIf flag = 2 Then
XX = midX
Else
XX = maxX
End If
e = 0.0000001 'epsilon
MinY = minY0
midY = midY0
maxY = maxY0
Do
LPmin = 0
LPmid = 0
LPmax = 0
For I = 1 To n
LPmin = LPmin - (Y(I) - MinY) / Sqr((X(I) - XX) ^ 2 + (Y(I) - MinY) ^ 2)
LPmid = LPmid - (Y(I) - midY) / Sqr((X(I) - XX) ^ 2 + (Y(I) - midY) ^ 2)
LPmax = LPmax - (Y(I) - maxY) / Sqr((X(I) - XX) ^ 2 + (Y(I) - maxY) ^ 2)
Next I
If Abs(LPmid) < e Then
Exit Do
ElseIf LPmid * LPmin < 0 Then
maxY = midY
Else
MinY = midY
End If
midY = (MinY + maxY) * 0.5
Loop
Worksheets("test").Cells(j, 7).Value2 = XX
Worksheets("test").Cells(j, 8).Value2 = midY
Worksheets("test").Cells(j, 11).Value2 = LPmid
L = 0
LPmin = 0
LPmid = 0
LPmax = 0
For I = 1 To n
L = L + Sqr((X(I) - minX) ^ 2 + (Y(I) - midY) ^ 2)
LPmin = LPmin - (X(I) - minX) / Sqr((X(I) - minX) ^ 2 + (Y(I) - midY) ^ 2)
LPmid = LPmid - (X(I) - midX) / Sqr((X(I) - midX) ^ 2 + (Y(I) - midY) ^ 2)
LPmax = LPmax - (X(I) - maxX) / Sqr((X(I) - maxX) ^ 2 + (Y(I) - midY) ^ 2)
Next I
If Abs(LPmid) < e Then
flag = -1
ElseIf LPmid * LPmin < 0 Then
maxX = midX
Else
minX = midX
End If
midX = (minX + maxX) * 0.5

Worksheets("test").Cells(j, 9).Value2 = L
Worksheets("test").Cells(j, 10).Value2 = LPmid
End Sub

Sub 最短距離2()
e = 0.0000001 'epsilon

'muuming2001様の複数のポイントがA,B,C,D,E　の5点でしたから N=5　としました。
n = 5
ReDim X(n) As Variant
ReDim Y(n) As Variant
'　EXCELシートからx,yのデータを読み取ります。
For I = 1 To n
X(I) = Worksheets("test").Cells(I + 1, 2).Value2
Y(I) = Worksheets("test").Cells(I + 1, 3).Value2
Next I

'x,yの最小値を決定します。
minX0 = X(1)
maxX0 = X(1)
For I = 1 To n
If minX0 > X(I) Then minX0 = X(I)
If maxX0 < X(I) Then maxX0 = X(I)
Next I
minX0 = minX0 * 1.01
maxX0 = maxX0 * 1.01
midX0 = (minX0 + maxX0) * 0.5

minY0 = Y(1)
maxY0 = Y(1)
For I = 1 To n
If minY0 > Y(I) Then minY0 = Y(I)
If maxY0 < Y(I) Then maxY0 = Y(I)
Next I
minY0 = minY0 * 1.01
maxY0 = maxY0 * 1.01
midY0 = (minY0 + maxY0) * 0.5

j = 1
Worksheets("test").Cells(j, 7).Value2 = "X"
Worksheets("test").Cells(j, 8).Value2 = "Y"
Worksheets("test").Cells(j, 9).Value2 = "L"
Worksheets("test").Cells(j, 10).Value2 = "dL/dx"
Worksheets("test").Cells(j, 11).Value2 = "dL/dy"
j = j + 1
'''' minX
flag = 1
minX = minX0
midX = midX0
maxX = maxX0
Call bisection(j, n, X, Y, minX, midX, maxX, minY0, midY0, maxY0, flag)
j = j + 1
'''' maxX
flag = 3
minX = minX0
midX = midX0
maxX = maxX0
Call bisection(j, n, X, Y, minX, midX, maxX, minY0, midY0, maxY0, flag)
j = j + 1
minX = minX0
midX = midX0
maxX = maxX0
Do
'''' midX
flag = 2
Call bisection(j, n, X, Y, minX, midX, maxX, minY0, midY0, maxY0, flag)
j = j + 1
If flag = -1 Then Exit Do
Loop
End Sub

これだと
x=4.12071877241135
y=6.56896227404475
でL=22.5247703740236となります。

回答No.7

f272

この問題は#3さんがいうように結構面倒でして，現実的には数値解法で求めるのが手っ取り早いでしょう。
その定式化としては，求める点を(x,y)とすれば
√((x-1)^2+(y-15)^2)
+√((x-4)^2+(y-2)^2)
+√((x-3)^2+(y-7)^2)
+√((x-9)^2+(y-10)^2)
+√((x-5)^2+(y-5)^2)
を最小化すればよいということになります。
真面目にといてみると
x=4.12071817415932
y=6.56896246705083
で最小値22.52477037になります。#4さんのやっているのは2分法もどきですから，いくら精度を高くしても真の解には近づきません。

お礼

みなさんありがとうございます。

コピペお礼で申し訳ありません。

かなり難しい問題なのがわかりました。
時間をかけて勉強してみたいと思います。

回答No.6

178-tall

参照 URL
　　↓
フェルマ点証明の解析的方法

参考URL：

http://aozoragakuen.sakura.ne.jp/taiwa/taiwaNch03/node7.html

回答No.5

noname#222520

ANo.1の回答者です。
ご指摘により、自分の考え方の誤りに気付きました。
異なる結果（距離の合計の最小値）からは、小数第一位を四捨五入すると等しくなるので、誤差の範囲程度にしか思っていませんでしたが、次の例から誤りは明白であると分かりました。

4^2+5^2+6^2+7^2+8^2(=190)＜1^2+2^2+3^2+9^2+10^2(=195)ですが、
1+2+3+9+10(=25)＜4+5+6+7+8(=30)であり、
それぞれの数を2乗したものの和と2乗しないままの和では、大小関係が逆転することもあります。

お礼

みなさんありがとうございます。

コピペお礼で申し訳ありません。

かなり難しい問題なのがわかりました。
時間をかけて勉強してみたいと思います。

回答No.4

noname#252159

No.2で数値解析したフローチャートとブログラムを提示します。
×が最適値に近づくと、次にyが最適値に近づくというループです。
ループの番人として、エプシロンを設けています。
下のVBA はExcel で作りました。動作確認済みですが、
　シート名をtest とすること
　　データ配列は変数xをA列に、変数yをB列において下さい。
　　　データ数Nはプログラムの部分で変更して下さい。

　プログラムはエレガントさに欠けます。いつか、より美しいプログラムを作りたいと思っています。
一-￣￣￣￣￣￣^-^

Sub 最短距離()
Dim X(100) As Variant
Dim Y(100) As Variant
'muuming2001様の複数のポイントがA,B,C,D,E　の5点でしたから N=5　としました。
N = 5
'　EXCELシートからx,yのデータを読み取ります。
For I = 1 To N
X(I) = Worksheets("test").Cells(I + 1, 2).Value2
Y(I) = Worksheets("test").Cells(I + 1, 3).Value2
Next I

'x,yの最小値を決定します。
MinX = X(1)
MaxX = X(1)
For I = 1 To N
If MinX > X(I) Then
MinX = X(I)
End If
If MaxX < X(I) Then
MaxX = X(I)
End If
Next I
MidX = (MinX + MaxX) / 2

MinY = Y(1)
MaxY = Y(1)
For I = 1 To N
If MinY > Y(I) Then
MinY = Y(I)
End If
If MaxY < Y(I) Then
MaxY = Y(I)
End If
Next I
MidY = (MinY + MaxY) / 2

'エプシロンとして、距離の精度を設定します。
e = 0.0001 'epsilon
J = 0

Do ' 下のLoop While まで、ループ計算します。
Lmin = 0
Lmax = 0
Lmid = 0
'x座標を２分法で、考えるxの範囲を絞ります。
For I = 1 To N
Lmin = Lmin + Sqr((X(I) - MinX) ^ 2 + (Y(I) - MinY) ^ 2)
Lmax = Lmax + Sqr((X(I) - MaxX) ^ 2 + (Y(I) - MaxY) ^ 2)
Lmid = Lmid + Sqr((X(I) - MidX) ^ 2 + (Y(I) - MidY) ^ 2)
Next I
dfMin = Lmin - Lmid
dfMax = Lmax - Lmid
If dfMin > dfMax Then
MinX = MidX
ElseIf dfMin < dfMax Then
MaxX = MidX
End If
MidX = (MinX + MaxX) / 2

Lmin = 0
Lmax = 0
Lmid = 0
'y座標を２分法で、考えるxの範囲を絞ります。
For I = 1 To N
Lmin = Lmin + Sqr((X(I) - MinX) ^ 2 + (Y(I) - MinY) ^ 2)
Lmax = Lmax + Sqr((X(I) - MaxX) ^ 2 + (Y(I) - MaxY) ^ 2)
Lmid = Lmid + Sqr((X(I) - MidX) ^ 2 + (Y(I) - MidY) ^ 2)
Next I
dfMin = Lmin - Lmid
dfMax = Lmax - Lmid
If dfMin > dfMax Then
MinY = MidY
ElseIf dfMin < dfMax Then
MaxY = MidY
End If
MidY = (MinY + MaxY) / 2
'次のMidX と　MidYが求める座標をワークシートに表示
Worksheets("test").Cells(J + 1, 5).Value2 = MidX
Worksheets("test").Cells(J + 1, 6).Value2 = MidY
J = J + 1
Loop While (dfMin > e Or dfMax > e)
'最小となったときの長さをワークシートに表示
Worksheets("test").Cells(J + 2, 6).Value2 = Lmid
End Sub

回答No.3

178-tall

3 ポイントないし凸四角形頂点なら、

　　　　　　　　↓

(参考URL) 【１】フェルマーの問題

それ以上は難問か。(数値解法など？)

参考URL：

http://www.geocities.jp/ikuro_kotaro/koramu/316_s.htm

お礼

みなさんありがとうございます。

コピペお礼で申し訳ありません。

かなり難しい問題なのがわかりました。
時間をかけて勉強してみたいと思います。

回答No.2

noname#252159

残念ながら、revenge-goさんの考え方では、最小にならないようです。
それよりも、小さい座標が存在します。

　muuming2001さんの例でいけば、
revenge-goさんの平方完成式は
5(x - 22/5)^2 +5(y-39/5)^2 + c　　　cは定数
座標(22/5, 39/5)　が求める点でその時の、
５点への距離の和は 23.35122
となります。
２乗したままの和で平方根の和に転化することはできないのです。

　しかし、２分法による計算では
求める座標は近似値で
x= 4.5625、
y= 6.164063
で５点への距離の和は22.58402
となります。

Program での2分法では如何

お礼

みなさんありがとうございます。

コピペお礼で申し訳ありません。

かなり難しい問題なのがわかりました。
時間をかけて勉強してみたいと思います。