已知橢圓5個點，用數學軟體求橢圓周長？

01-16

已知橢圓五個點坐標如下
求橢圓的周長

http://mathworld.wolfram.com/ConicSection.html

Mathematica代碼

pts = {{5.764, 0.648}, {6.286, 1.202}, {6.759, 1.823}, {7.168, 2.526}, {7.48, 3.36}}; eq = Function[{x,y}, {x^2, x y, y^2, x, y, 1}] @@@ Append[pts, {x,y}]//Det//Factor ContourPlot[eq == 0, {x, -2, 8}, {y, -2, 8}, Epilog -&> Point@pts] ArcLength@ImplicitRegion[eq == 0, {x, y}] // InputForm

======================================================

不用ArcLength的方法：

pts=Rationalize@{{5.764,0.648},{6.286,1.202},{6.759,1.823},{7.168,2.526},{7.48,3.36}}; eq=Function[{x,y},{x^2,x y,y^2,x,y,1}]@@@Append[pts,{x,y}]//Det//Factor//Last {a,b}=Module[{a,b,c,d,e,f,g},{f,{d,e},{{a,b},{g,c}}}=CoefficientArrays[eq,{x,y}]; {b,d,e}={b,d,e}/2; Sqrt[(2 Det[{{a,b,d},{d,e,f},{b,c,e}}])/(Det[{{a,b},{b,c}}] (a+c+Sqrt[(a-c)^2+4 b^2] {-1,1}))]]; ContourPlot[{eq==0,x^2/a^2+y^2/b^2==1},{x,-5,8},{y,-5,8},Epilog-&>Point@pts] 4 a EllipticE[1-b^2/a^2] N[%,30]

還可以通過解方程求a, b

sys=Function[{x,y},{1/aa,1/bb}.(RotationMatrix[t].{x-x0,y-y0})^2-1]@@@pts sol=NSolve[sys==0,{aa,bb,x0,y0,t}]

更新：沒想到真的被「個人任務」了，我就說下不是我舉報的啊

--------------------------------------------------

這個應該算是數學建模中比較常見的問題，本來以為這題會被「個人任務」的，沒想到反而關注的人還挺多

yellow的想法不錯，可以求出一般方程的係數：

p = [5.764 0.648;6.286 1.202;6.759 1.823;7.168 2.526;7.48 3.36]; f = @(x,y)[x.^2 x.*y y.^2 x y x==x]; a = null(f(p(:,1),p(:,2))); plot(p(:,1),p(:,2),"o"), hold on h = ezplot(@(x,y)f(x(:),y(:))*a,[-2 8]);

不過只從一般方程的係數並不好得到曲線周長，而且MATLAB又沒有類似於Mathematica的ArcLength這樣內置函數，所以還是考慮求出其對應的標準型來求橢圓積分比較方便。

首先通過x和y項的係數求出橢圓心坐標，並將數據平移至原點：

c = [a(4) a(5)]/[2*a(1) a(2);a(2) 2*a(3)]"; plot(c(1),c(2),"r+"); p = bsxfun(@plus,p,c);

之後根據旋轉後的標準方程求出a,b：

f = @(a,b,t)(p*[cos(t);sin(t)]).^2/a^2+(p*[sin(t);-cos(t)]).^2/b^2-1; B = fminsearch(@(p)norm(f(p(1),p(2),p(3))),rand(3,1),optimset("TolX",1e-9));

之後可以用橢圓積分算出周長:

&>&> l = 4*B(1)*integral(@(r)sqrt(1-(1-B(2)^2/B(1)^2)*sin(r).^2),0,pi/2) l = 21.8344706519725

圖就不放了，和上邊幾個都差不多，而且可以自己運行得到

上述過程中提前平移到原點的處理主要是為了後邊計算更穩定，其實可以把橢圓心也當做參數來一起求，完整的程序是：

p = [5.764 0.648;6.286 1.202;6.759 1.823;7.168 2.526;7.48 3.36]; B = lsqnonlin(@(B)(bsxfun(@minus,p,B(1:2))*[cos(B(5));sin(B(5))]).^2/B(3)+... (bsxfun(@minus,p,B(1:2))*[sin(B(5));-cos(B(5))]).^2/B(4)-1,ones(1,5)); l = 4*sqrt(B(3))*integral(@(r)sqrt(1-(1-B(4)/B(3))*sin(r).^2),0,pi/2)

得到：

l = 21.8344705313394

不過這樣穩定會有所下降，偶爾並不能首次就得到正確結果，如果只計算一次，可以人為判斷優化函數的其餘輸出值（例如exitflag）的值，如果需要多次判斷可以寫一個循環

另外還有一種情況就是當精度要求不高時可以直接根據作圖用的曲線求出長度：

&>&> l = sum(sqrt(sum(diff(h.ContourMatrix(:,2:end),[],2).^2,1))) l = 21.8339855020675

這裡的h是第一段代碼中ezplot所返回的句柄

另外，根據向量場相關支持可知道，對於二維區域：

$D = left{(x,y):f(x,y)<=0 ight}$

其邊界為：

$C=partial D = left{(x,y):f(x,y)=0 ight}$

則有：

$iint_D ablacdot{f v} ds = oint_C{f v cdot n} dell$

為求C曲線長，可以取：

${f v = n} = frac {( abla f)}{|| abla f||}$

使得：

$iint_D ablacdot frac {( abla f)}{|| abla f||} ds = oint_C dell$

所以假設橢圓方程是f(x,y)=0, 顯然在橢圓內部有f(x,y)&<0, 於是可以通過符號計算計算上述量，並且最後做在f(x,y)&<=0區域上的二重積分即可得到周長

不過其中過程需要符號計算，如果MATLAB符號計算工具箱的話可以嘗試一下，但是性能可能會比較差或者最後一步積分的數值計算精度有問題，所以還是比較推薦上邊的純數值方法，符號計算這種事情還是交給Mathematica或者Maple來做比較方便

瀉藥，這個問題我嘗試了一下，但一直得不到理想的答案，這裡就當為MATLAB用戶拋磚引玉了。

首先，我們來看對於2-D平面內的一個橢圓，其一般方程如下：

$f(x, y)=Ax^2+Bxy+Cy^2+Dx+Ey+Fequiv 0$

為何一定有如上形式表達，其實可以這樣理解，任意一個平面內的橢圓都是由一個標準橢圓： $frac{x^2}{a^2}+frac{y^2}{b^2}=1$ 經過平移和旋轉得到，相當於在標準方程上對 $(x, y)$ 進行線性變換，其具體推導後面會提及，暫且按下不表。

現在根據已知的5個數據點 $(x_1, y_1), (x_2, y_2), cdots, (x_5, y_5)$ ，將其代入該方程應當滿足如下齊次線性方程組：

$left( egin{array}{cccccc} x_1^2 x_1y_1 y_1^2 x_1 y_1 1\ x_2^2 x_2y_2 y_2^2 x_2 y_2 1\ vdots vdots vdots vdots vdots vdots\ x_5^2 x_5y_5 y_5^2 x_5 y_5 1\ end{array} ight) left( egin{array}{c} A\ B\ C\ D\ E\ F end{array} ight) = extbf{0}$

為方便表述，進行如下簡記：

$A = left( egin{array}{cccccc} x_1^2 x_1y_1 y_1^2 x_1 y_1 1\ x_2^2 x_2y_2 y_2^2 x_2 y_2 1\ vdots vdots vdots vdots vdots vdots\ x_5^2 x_5y_5 y_5^2 x_5 y_5 1\ end{array} ight), x = left( egin{array}{c} A\ B\ C\ D\ E\ F end{array} ight)$

不失一般性，假設 $lVert x Vert = 1$ (即將橢圓方程係數向量 $x$ 進行單位化)，於是以上方程變為一個帶約束的欠定齊次方程組最小二乘問題：

$min_x~lVert Ax Vert_2^2~~s.t.~lVert x Vert=1$

這個問題的求解可以藉助SVD（Singular Value Decomposition），我們對 $A$ 進行相應分解得到 $A = USigma V^*$ ，其中 $U, V$ 為酉矩陣（其各列分別為 $A^*A, A^*A$ 的特徵向量）， $Sigma$ 為一個對角矩陣，其對角線元素從大到小排列其奇異值 $sigma_igeq 0$ ，那麼當且僅當 $x$ 為 $V$ 的最後一列時（最小奇異值對應的右奇異向量） $lVert Ax Vert_2^2$ 達到最小值 $sigma_{min}^2$ .按照這個思路，很快能得到橢圓的軌跡方程：

$0.1690x^2-0.1892xy+0.1910y^2-0.8996x-0.1191y+0.2751equiv0$

%% x = [5.764, 6.286, 6.759, 7.168, 7.480]"; y = [0.648, 1.202, 1.823, 2.526, 3.360]";


M = [x.^2, x.*y, y.^2, x, y, ones(size(x))];

[U S V] = svd(M);
scatter(x, y, "SizeData", 40, "MarkerFaceColor", "b", "MarkerFaceAlpha", .75)

hold on

h = ezplot(@(x, y) V(1, 6)*x.^2+V(2, 6)*x.*y+V(3, 6)*y.^2+V(4, 6)*x+V(5, 6)*y+V(6,6), ...

       [-2, 10, -2, 6]);

h.set("LineWidth", 1.5)
xlabel("$x$", "interp", "latex", "FontSize", 13)

ylabel("$y$", "interp", "latex", "FontSize", 13)

title("A fitting ellipse using SVD and LLS", "FontSize", 14)

grid on

%% saveas(gcf, "ellipse.png")

但這個結果還遠沒達到我預想的地步----我希望通過對一個標準橢圓進行平移和旋轉得到該目標橢圓，這樣我能獲得更多額外的信息，例如橢圓中心點坐標 $(x_0, y_0)$ ，半長短軸 $a, b$ 以及對應的旋轉角度 $heta$ .

一旦這些信息被獲得，那麼計算橢圓周長只需要通過 $a,b$ 來計算即可。雖然橢圓的周長不存在一個初等的公式表達，但對標準橢圓周長計算的研究結果（橢圓積分）要比通過對閉合曲線直接進行分割計算周長來得可靠。

最後說說我最開始的設想吧，我們假設該橢圓是由一個標準橢圓（焦點在 $x$ 軸），長半軸和短半軸長分別為 $a, b$ ，將其逆時針旋轉 $heta$ 後，再將中心平移到 $(x_0, y_0)$ 處。那麼標準橢圓上一點 $(x, y)$ 經如此變幻後其坐標變為 $(x$ ，不難得到如下等式：

$left( egin{array}{c} x$

於是很快可以得到 $(x$ 的軌跡方程：

$b^2Big(cos( heta)(x$

PS.根據上述方程就能佐證之前我們對一般橢圓都具有二次形式的斷言。

代入五個點求解五個參數 $(a^2, b^2, x_0, y_0, heta)$ ，其中 $a > b$ ，直接通過fsolve/solve時無法得到解的，需要藉助一些優化技巧，具體求解過程可以參照Falccm的答案--已知橢圓5個點，用數學軟體求橢圓周長？ - Falccm 的回答。

-------------------------------------- &>&>&> 2016/02/20最後更新 &<&<&<--------------------------------------

根據 @曹洪洋的建議（在此表示感謝），通過將 $cos( heta), sin( heta)$ 轉為自由變數，同時加入等式約束 $cos( heta)^2+sin( heta)^2=1$ ，這樣相當於6個方程求解6個變數。經嘗試利用solve就能求解。

求解出 $a, b$ 後，根據橢圓標準方程計算其周長，這裡採用參數方程形式進行積分：

$4int_0^{frac{pi}{2}}sqrt{a^2sin( heta)^2+b^2cos( heta)^2}d heta$

%% x = [5.764, 6.286, 6.759, 7.168, 7.480]"; y = [0.648, 1.202, 1.823, 2.526, 3.360]";


syms a b x0 y0 c s

f = arrayfun(@(x, y) b*(c*(x-x0)+s*(y-y0))^2+a*(-s*(x-x0)+c*(y-y0))^2-a*b, ...

             x, y, "UniformOutput", false);
rs = solve([f{:} c^2+s^2-1], [a, b, x0, y0, c, s]);

a = double(rs.a); a = a(a~=0); % non-zero solutions

[a2 b2] = deal(max(a), min(a));

perimeter = 4*integral(@(t)sqrt(sin(t).^2*a2+cos(t).^2*b2), 0, pi/2)
perimeter =

21.8345

參考資料：

[1] How to plot an Ellipse

[2] 應用奇異值分解求解最小二乘法問題

[3] 線性最小二乘問題

不用版本10函數的且更粗暴一些的做法：

eqn = 0 == 1 + {x^2, x y, y^2, x, y}.a /@ Range[5] expr = y /. Solve[eqn, y] /. Solve[({x, y} [Function] #) @eqn @@@ pts][[1]] delta = Cases[expr, Sqrt[_], ∞][[1]]


Plot[expr, {x, -10, 10}]

2 NIntegrate[Boole[delta &> 0] Sqrt[1 + D[expr[[1]], x]^2], {x, -10, 10}]

21.8345