問題描述
在平面內有n（n>=3）個點N1(x1,y1),N2(x2,y2),...,Nn(xn,yn)，現求一點P(x,y)，使得P到各點直線距離之和最小。
算法分析
當n=3時，這是著名的三角形費馬點問題，網上有詳細介紹和證明。
然而，那些平面幾何證明看似巧妙，但真正涉及到了n個點的時候，就只能呵呵了，還是得用解析法來想辦法。
目標函數為：
我們需要求它的最小值。
分別對x和y求偏導數：
fx(x,y) =
=0
fy(x,y) =
=0
當兩偏導數同為0的時候，是此二元函數的駐點。
這里，若再求一次偏導數（二階偏導數）：
發現其恆大於0，即原函數是個凸函數，一階偏導數為0的點就是它的最小值點。
那么上面這倆方程（一階偏導數為0）怎么解呢？
這里就要隆重推出：迭代法。
說句題外話，我們用通俗的說法來辨析幾個詞：循環、迭代、遞歸、遍歷。我們常常把它們混用，但它們其實是互不相同的！
循環Repeating：反復執行一段程序代碼。例如：while循環、for循環...
迭代Iteration：每一次的運算結果會成為下一次運算的初始值，經常用此方法以不斷逼近目標結果。例如：牛頓迭代法、二分法、斐波那契數列...
遞歸Recursion：（函數等）自己調用自己。例如：漢諾塔問題、斐波那契數列...
遍歷Traversal：訪問一個樹的所有結點，每個只訪問一次。例如：廣度優先搜索、深度優先搜索...
可以看出，有些問題可能會同時涉及到這四個中的多個（其實也不難理解），這大概就是我們常常混淆這四個詞的一大原因吧。比如，實現后三者的程序，寫代碼時基本都跑不掉循環結構，很多人從此就開始混淆它們了...而其實，很多時候它們的思想也是相通的，因此造成了一個問題既可以用迭代也可以用遞歸的現象，這很正常。
回到本文的問題。我們嘗試解一個方程f(x)=0，如果能找到一個g(x)=f(x)+x，從而將原方程轉化為x=g(x).通過不斷迭代：g(g(g(g(g(x)))))，逼近解x。
這其實就是高中數學（一般在競賽中出現，或是十多年前的很難的高考數學的數列題里出現）里的“不動點”問題。
一個經典的例子就是利用cos(cos(cos(0)))解方程cos(x)=x。利用作圖能清晰的看出：
當然了，此種方法對g(x)的圖形和x的初始值是有要求的，不然有可能不但不逼近，反而跑的越來越遠了。
例如，對於方程3^x-7=x，我們直觀的感覺它有兩個解，運行如下代碼，輕松得到了一個解-6.9995。
import math
m=0
for i in range(10):
m=3**m-7
print(m)
可另外一個解呢？此種方法就很難求了，因為一旦迭代就離它越來越遠了。我們嘗試更改初始值如m=7來求那個正數的解：
根本執行不了。其實在草稿紙上畫畫圖就能發現，想得到那個解，通過我們這種方法根本不能收斂。
不扯遠了，再次回到本題，將x和y可以表示為：
初始值可設為那n個點的重心，其實，最終答案一般都不會離重心太遠。
至於為什么本題就可以收斂？由二階導數>0恆成立知，一階導數單調，故上述方程僅有一個根。
我們可以作圖看看右邊這個奇怪的函數和函數f(x)=x的交點，運行以下代碼：
import math
import random
import matplotlib.pyplot as plt
a1=[];a2=[]
for i in range(6):
a1.append(random.random())
a2.append(random.random())
for j in range(1000):
x=j/1000;y=j/1000
xfenzi=0;xfenmu=0;yfenzi=0;yfenmu=0
for i in range(6):
g=math.sqrt((x-a1[i])**2+(y-a2[i])**2)
xfenzi=xfenzi+a1[i]/g
xfenmu=xfenmu+1/g
yfenzi=yfenzi+a2[i]/g
yfenmu=yfenmu+1/g
xn=xfenzi/xfenmu
yn=yfenzi/yfenmu
plt.scatter(x,xn,color='b',s=1)
plt.scatter(x,x,color='g',s=1)
plt.show()
由上述，顯然，不論怎么運行，仍是只有一個根。而且很明顯這個圖形是符合我們描述的那個迭代的，即：通過迭代能逐漸逼近那個交點。
因此，我們可以用迭代法實現找出平面上到n個點距離之和的最小值的費馬點了。
以下是Python實現代碼：
import math
import random
import matplotlib.pyplot as plt
n=int(input('請輸入n：'))
a1=[];a2=[]
for i in range(n):
a1.append(random.random())
a2.append(random.random())
plt.scatter(a1,a2,color='r')
x=sum(a1)/n;y=sum(a2)/n
while True:
xfenzi=0;xfenmu=0;yfenzi=0;yfenmu=0
for i in range(n):
g=math.sqrt((x-a1[i])**2+(y-a2[i])**2)
xfenzi=xfenzi+a1[i]/g
xfenmu=xfenmu+1/g
yfenzi=yfenzi+a2[i]/g
yfenmu=yfenmu+1/g
xn=xfenzi/xfenmu
yn=yfenzi/yfenmu
if abs(xn-x)<0.01 and abs(yn-y)<0.01:
break
else:
x=xn
y=yn
plt.scatter(x,y,color='b')
plt.show()
運行效果：
