おはようございます。週明けから、保育園から呼び出しをくらった熱の高めの長男と実家で過ごして、戻ってきた@kjunichiです。

背景

たまたまTLで以下を見かけた。

http://togetter.com/li/1036604

1/7がつくる楕円の方程式

1/7=0.142857142857142857...

と1,4,2,8,5,7が繰り返し出てくる。ここから点(1, 4), (4, 2), (2, 8), (8, 5), (5, 7) の組み合わせ選ぶと、これらの点を通る楕円が描けるという。

楕円の方程式というか、二次曲線の方程式の一般形は、

a*x^2+b*x*y + c*y^2 + d*x + e*y +f = 0

で表せるとのこと。今回描画したい楕円もa,b,c,d,e,fが決まればJuliaでプロットできる。真のJulia使いなら、この５つのパラメータを求めるところからなのだろうが、ヘタレなので、 Mathematicaで有名なWolframリサーチ社の力を借りる。

http://mathworld.wolfram.com/One-SeventhEllipse.html

その結果、以下の数式を得ることが出来る。

19* x^2+36*y*x+41*y^2-333*x-531*y+1638

Juliaでやる

プロットする数式が求まったので、以降はJuliaで進める。

必要なもの

陰関数をプロットするライブラリ

f(x,y)=0を満たすｐ点をプロットする便利なあライブラリがった。

Pkg.add("Plots")
Pkg.add("ImplicitEquations")

楕円をプロットする

using Plots
pyplot()
using ImplicitEquations
f(x,y) =  19* x^2+36*y*x+41*y^2-333*x-531*y+1638
plot(f == 0, (-2,10),(0,10))

残念ながら、プロットするx,yの範囲を指定しないといい感じの楕円の。グラフにならない。

f:id:kjw_junichi:20161029075206p:plain

Gadflyでもやってみた

変数の範囲を指定するのであれば、Gadflyでも出来るだろうという事でやってみた。

using Gadfly
f(x,y) =  19* x^2+36*y*x+41*y^2-333*x-531*y+1638
plot(z=(x,y)->f(x,y),
     x=linspace(-1,10,100), y=linspace(0,10,100), Geom.contour(levels=[0]))

f:id:kjw_junichi:20161029075234p:plain

やっぱりJuliaで一般形からa,b,c,d,e,fの各係数を求めてみた

fは定数なので、1と、仮定して、連立方程式から行列を作る。

A = [
1//1 4//1 16//1 1//1 4//1;
16//1 8//1 4//1 4//1 2//1;
4//1 16//1 64//1 2//1 8//1;
64//1 40//1 25//1 8//1 5//1;
25//1 35//1 49//1 5//1 7//1;
]

//1とすることで、分数表示で以降の計算を。行ってくれる。（さらっと書いたが、これに気が付くのに半日以上かかった。。）

B=inv(A)
C=B*[1;1;1;1;1]

5-element Array{Rational{Int64},1}:
 -19//1638
  -2//91  
 -41//1638
  37//182 
  59//182

分母に1638が出現するf=1638と

C*1638

5-element Array{Rational{Int64},1}:
   19//1
   36//1
   41//1
 -333//1
 -531//1

学んだこと

任意の5点指定すれば、それらを通る楕円が１つに決まる。(必ずしも楕円になるとは限らない)
Julia言語でも陰関数をプロットするライブラリがあった。
陰関数のプロットに関して論文が出ていることを知った。（http://www.dgp.toronto.edu/people/mooncake/papers/SIGGRAPH2001_Tupper.pdf）
行列の成分をN//1とすると、逆行列が分数の形で取得できる.