mikeo_410
  1. 2.地球
    1. 1.地球楕円体上の点
    2. 2.地球儀と楕円体の描画
    3. 3.3Dの直線と平面
    4. 4.卯酉線
    5. 5.楕円の標準形
    6. 6.5点で決まる楕円
    7. 7.曲率半径と平均曲率
    8. A1.地図、地球儀(Rスクリプト、データ)
    9. A2. 付録6.計算式集の公式

楕円の標準形

  1. 楕円を描く
  2. 回転、平行移動した楕円の式
  3. 楕円の諸元から一般形を求める
  4. 一般形から楕円の諸元を求める
  5. 一般形から楕円の諸元を求める例
  1. 最初の例
  2. 90°の傾き
  3. 縦長の標準形
  4. 回転角の計算
  5. ゼロ判定と精度
  6. 定数項Fの役割
  7. 平方完成
  8. 回転と係数の関係
  1. 標準形との点の対応
  1. 標準形からグラフを描く
  2. 一部を描く
  3. 座標指定で標準形の弧を描く
  4. 標準形から楕円体表面へ
  1. 行列の固有値との関係
  1. 固有値
  2. 固有ベクトル
  1. 固有値による標準形の計算
  1. 係数A’、C’
  2. 係数D’、E’
  1. eigen関数による計算例
  2. 固有値による標準形の計算の図示
  3. 固有ベクトルによる座標変換
  4. 行列による楕円の描画
  5. 解法の逆順での楕円の描画
  6. 平方完成
  7. 都合の良い固有値、固有ベクトル
  1. 固有値、固有ベクトルの計算
  2. 都合の良い計算
  3. eigen関数の結果の組み替え

1. 楕円を描く

実際に楕円を描くには媒介変数表示を使用していると思います。

  1. a <- 3; b <- 2
  2. θ <- seq(0,2*pi,length.out=361)
  3. x <- a*cos(θ)
  4. y <- b*sin(θ)
  5. plot(x,y,t="l",xlim=c(-3,3),ylim=c(-3,3))
  6. abline(h=c(0),v=c(0),col="gray40")

これは左図のように、角度 θ に対して、点Pを描くことです。

 x軸半径 a 、y軸半径 bの楕円を考えます。点Qは半径aの補助円状にあります。線分OQと、x軸の成す角がθです。このとき、楕円上の点P(a・cosθ、b・sinθ) を描きます。

原点O と点P を結んだ直線が、x軸となす角を γ とすると、   なので、

    

    

 後者を前者で割って、

   

 したがって、

     (1.1)

     (1.2)

 この関係は、一部を描く場合に役立ちます。

下図の、楕円の周上の2点P1(x1、y1)、P2(x2、y2) の間の弧を描くには、その媒介変数 θ1、θ2 を求めます。まず、中心角は、それぞれ、 です。

したがって、θ1、θ2 は、 です。

  1. plot(0,0,pch=16,xlim=c(-3,3),ylim=c(-3,3))
  2. abline(h=0,v=0,col="gray60")
  3. a <- 3; b <- 2
  4. t <- seq(0,2*pi,length.out=361)
  5. lines(a*cos(t),b*sin(t),col="cyan")
  6. x1 <- 1; y1 <- (b*sqrt(a^2-x1^2))/a
  7. x2 <- 2; y2 <- (b*sqrt(a^2-x2^2))/a
  8. abline(a=0,b=y1/x1,col="gray80")
  9. abline(a=0,b=y2/x2,col="gray80")
  10. abline(v=c(x1,x2),h=c(y1,y2),col="gray80")
  11. θ1 <- atan(a*y1/(b*x1)); θ2 <- atan(a*y2/(b*x2))
  12. t <- seq(θ1,θ2,by=(θ2-θ1)/20)
  13. lines(a*cos(t),b*sin(t),lwd=5,col="blue")

x軸半径を a、y軸半径を b として、楕円は、

                        (1.3)

と表され、θは媒介変数です。もし、a=r、b=r なら、

                        (1.4)

半径 r の円を描きます。

また、

                                        (1.5)

は、半径1の円を描きます。円周上の点のX座標とY座標の2乗の和は、原点からの距離の2乗です。

                    (1.6)

式(1.3) を、以下のように書き換えて、

式(1.6) に当て嵌めれば、楕円の式になります。

                                                (1.7)

円や楕円の基本は、

                                        (1.8)

で、 や、 のように考えられます。

また、式に定数を乗じたり、割ったりしても、解やグラフは変わりません。

 半径 r の円の方程式は、 や、 と表されます。楕円の式も、分数を避けて、

                        (1.9)

と表せます。

式(1.7) は、原点を中心とした楕円を表しますが、原点を (m、n) に移動する平行移動をした式は、

と置いた、式(1.8) だと考えることができます。

点 (x、y) を原点の周りに反時計回りにβ回転すると、(x・cosβ+y・sinβ、 -x・sinβ+y・cosβ) となります。

β回転した楕円の式は、

        

と置いた、式(1.8) だと考えることができます。

2. 回転、平行移動した楕円の式

x軸半径 a、y軸半径 b、の楕円を反時計回りに γ 回転し、楕円の中心を原点から(m、n)に移動した楕円の式を求めます。

回転していない楕円を、媒介変数θ を使って表すと、

回転後の楕円は、

さらに、x軸方向に m 、y軸方向に n 平行移動します。

となります。

これは、以下のように考えることができます。

1)単位円の座標列に、2)伸縮行列を乗じ、3)回転行列で回転して、4)平行移動します。

この媒介変数の部分を、

と置き換えると、

           (2.1a)

            (2.1b)

X、Y を求めると、

                  (2.2a)

                (2.2b)

となります。

また、式(2.1a)(2.1a) は、

また、式(2.2a)(2.2a) は、

と表せます。

左図は、

の楕円の例です。

  1. a <- 4; b <- 2; γ <- pi/6;
  2. m <- sqrt(3); n <- 2
  3. θ <- seq(0,2*pi,length.out=361)
  4. X <- cos(θ); Y <- sin(θ)
  5. xy <- rbind(c(cos(γ),-sin(γ)),
  6.             c(sin(γ),cos(γ))) %*%
  7.       rbind(c(a,0),c(0,b)) %*%
  8.       rbind(X,Y) + c(m,n)
  9. plot(xy[1,],xy[2,],t="l",col="blue")
  10. abline(v=c(0,sqrt(3)),h=c(0,2))
  11. abline(a=2-m*tan(γ),b=tan(γ))
  12. abline(a=2+m*(1/tan(γ)),b=-1/tan(γ))

スクリプトは、式(2.1a)(2.1a) の計算で、傾いた楕円を描きます。

また、楕円の周の座標の列に対して、式(2.2a)(2.2a) の演算を行うと単位円の座標に戻ります。

  1. XY <- rbind(c(1/a,0),c(0,1/b)) %*%
  2.       rbind(c(cos(γ),sin(γ)),
  3.             c(-sin(γ),cos(γ))) %*%
  4.       (xy - c(m,n))
  5. plot(XY[1,],XY[2,],t="l")
  6. round(range(XY[1,]-X),5)
  7. round(range(XY[2,]-Y),5)

下段の図のように、式(2.1a)(2.1a) は、(x、y)=(1、0) を楕円の周の赤丸(3√3、4) に変換します。

式(2.2a)(2.2a) の変換は、楕円の周の赤丸(x、y)=(3√3、4) を、(X、Y)=(1、0) に変換します。

式(2.2a)(2.2a) の変換は、xy座標の座標値で表された楕円の周の座標値を、XY座標で表した単位円の周の座標値に変換します。

また、xy座標で表した単位円の周の座標値に変換すると同じです。

半径1の円の式(1.8) は、以下のようになります。

                (2.3)

x軸半径 a、y軸半径 b、の楕円を反時計回りに γ 回転し、楕円の中心を原点から(m、n)に移動した楕円の周上の点(x、y) は、式(2.1) を満たします。

この式(2.1) に を乗じて整理します。

                                        (2.4)

2次曲線の一般形、

の係数は、        


(2.5)

楕円の式に実数を乗じても同じグラフになるので、同じグラフを描く楕円の式は無数にあります。

ここでは、 を乗じて、分数を避けた式になりました。

3. 楕円の諸元から一般形を求める

楕円の諸元を、

とします。

標準形の式で、

                        (3.1)

と表される楕円を、30°回転して、x軸方向に√3、y軸方向に2、平行移動した楕円の式を求めます。

式(2.5) によって、一般形の係数を求めることになります。

R だと、

  1. r$> GeneralFormsCoeffEllipse <- function(a,b,γ=0,m=0,n=0){
  2.       A <- a^2*sin(γ)^2+b^2*cos(γ)^2;
  3.       B <- 2*(b^2-a^2)*cos(γ)*sin(γ);
  4.       C <- b^2*sin(γ)^2+a^2*cos(γ)^2;
  5.       D <- -2*a^2*m*sin(γ)^2-2*(b^2-a^2)*n*cos(γ)*sin(γ)-2*b^2*m*cos(γ)^2;
  6.       E <- -2*b^2*n*sin(γ)^2-2*(b^2-a^2)*m*cos(γ)*sin(γ)-2*a^2*n*cos(γ)^2;
  7.       F <- (b^2*n^2+a^2*m^2)*sin(γ)^2+(2*b^2-2*a^2)*m*n*cos(γ)*sin(γ)+
  8.       (a^2*n^2+b^2*m^2)*cos(γ)^2-a^2*b^2;
  9.       list(A=A,B=B,C=C,D=D,E=E,F=F)
  10.     }
  11.     u <- GeneralFormsCoeffEllipse(4,2,pi/6,sqrt(3),2)
  12.     str(u)
  13. List of 6
  14.  $ A: num 7
  15.  $ B: num -10.4
  16.  $ C: num 13
  17.  $ D: num -3.46
  18.  $ E: num -34
  19.  $ F: num -27

maxima だと、

となって、楕円の一般形は、

                (3.2)

と求められました。

4. 一般形から楕円の諸元を求める

一般形の楕円の式が与えられたとき、その諸元を求めることを考えると、与式の係数には定数が乗じられていることを前提に考える必要があります。したがって、式(2.5) を、以下のように修正して解くことにします。U は任意の実数です。


(4.1)

 を求めたら、回転していない場合の係数を計算することになります。回転していないと言うのは楕円の軸と座標軸が平行だと言うことです。

γ=0 のとき、式(4.1) は、


(4.1a)

もし、γ=90° を採れば、式(4.1) は、


(4.1b)

与式は、2次曲線の一般形の、

                (4.2)

で、楕円の判定条件、を満たすものとします。

また、係数Aが負なら、式に -1 を乗じて、A>0 のケースだけを考えることにします。

式(4.1) の A、C の式から、

(4.3)

また、Bの式から、

                        (4.4)

式(4.1) の D、E に、式(4.3)、(4.4) を適用して、

連立方程式を解いて、

(4.5)

移動した楕円の中心 (cx、cy) は、元が (0,0) なので、(m,n) です。

与式の係数F は、式(2.3) の係数FにUを乗じたものなので、式(4.3)、(4.4) 、(4.5) を適用して、

 

 

                (4.6)

 なので、式(4.3) から、

                (4.7)

また、式(2.5) の係数B の式の両辺を2乗して、三角関数の部分を置き換えると、

        

                                        (4.8)

式(4.7)、(4.8) から、2組の  が得られます。

                (4.9)

                (4.10)

 とすると、

式(4.9) は、

                        (4.11)

式(4.10) は、

                        (4.12)

回転角度は、

                (4.13)

 と置き換えると、

                        (4.14)

式(4.11) からは、

                (4.15)

式(4.12) からは、

                (4.16)

2乗されていることから γ は、0から90°での値域で計算されることになります。楕円は180°回転すると重なるので、180°まで必要です。これは、係数Bの符号で判断できます。

式(4.6) から、

                        (4.17)

楕円の式、

が与えられた場合、式(4.11) 、あるいは、式(4.12) で  を求めることができます。

 なので、

                                (4.18)

                        (4.19)

と、楕円の半径 a、b を計算できます。

また、与式に対して、回転していない楕円(γ=0)の式は以下のように表せます。

平行移動した楕円は、

 を乗じて整理すると、

                                        (4.20)

                                        (4.21)

展開すると、

                (4.22)

5. 一般形から楕円の諸元を求める例

最初に係数Aが負なら、式に-1を乗じて、必ず、A>0 とすることにします。

5.1. 最初の例

例として、

を掲げます。係数は、

です。

式(4.11) から、

      (5.1.1)

となります。

式(4.15) から、

        

30°

※atan()の出力は、[-π、π] なので、0から180°まで連続した角度で表すには係数Bの符号を参照する必要があります。

ここでは式(4.11) を採ったので、式(4.4) の右辺の分母は負になります。係数Bも負なので、左辺 が正であることが分かります。

式(4.5)から、

式(4.17)、(4.18)、(4.19) から、

楕円の半径は、a=4、b=2 と分かります。

したがって、

は、X軸半径4、Y軸半径2の標準形の楕円を、反時計回りに π/6 回転し、X軸方向に√3、Y軸方向に2移動した楕円です。

あるいは、

は、中心が (√3、2) にある、30°傾いた、長半径が4、短半径が2 の楕円です。

5.2. 90°の傾き

2次曲線の一般形の式の係数から楕円の諸元を求めるスクリプトを載せます。そのためには1つ決めなければならないことがあります。標準形が縦長か横長かを決めて90°の傾きを使用するか、90°の傾きは2つの半径を入れ替えることにして90°の傾きは存在しないことにするか、と言うことです。

諸元から楕円の式を作る場合は、左図のように、

 

と、

 

の楕円は同じ式となり青い楕円になります。

 また、

 

と、

 

の楕円は同じ式となり赤い楕円になります。

 左図の2つの楕円は共に楕円の軸と座標軸が一致し、傾いていない楕円だと考えれば、90°の傾きの楕円はないことになります。違いは半径に現れます。

これに対して、青い楕円、a>b を標準形と定めれば、赤い楕円は90°傾いた青い楕円です。赤丸の半径は青丸と同じで、傾きが違う楕円です。

ここでは、後者を採ります。

したがって、ここに上げる StandardEllipse関数は、

   

の値を返します。

 左図は、 の条件で、 について、式(4.1)を計算したものです。

 係数Bがゼロになるのは、 γ=0 と γ=90° です。γ=0 のとき、係数Aは最小、係数Cは最大となります。γ=90° のとき、係数Aは最大、係数Cは最小となります。係数A,Cの最小値、最大値は同じ値で、 です。

StandardEllipse関数は、以下のようになります。

  1. StandardEllipse <- function(A,B,C,D,E,F){
  2.   if(A<0){A<- -A;B<- -B;C<- -C;D<- -D;E<- -E;F<- -F;}
  3.   if(B==0) {
  4.     if(A<C){ γ <- 0;     La <- C; Lb <- A;
  5.     } else { γ <- pi/2;  La <- A; Lb <- C; }
  6.   } else {
  7.     u <- sqrt(B^2+(C-A)^2); v<- C+A
  8.     La <- (u+v)/2; Lb <- (-u+v)/2
  9.     γ <- atan2(sqrt(u-(C-A)),-sign(B)*sqrt(u+(C-A)))
  10.   }
  11.   g <- 4*A*C-B^2
  12.   m <- (B*E-2*C*D)/g
  13.   n <- (B*D-2*A*E)/g
  14.   Ua2b2 <- (A*E^2-B*D*E+C*D^2)/g - F
  15.   list(a=sqrt(abs(Ua2b2/Lb)),b=sqrt(abs(Ua2b2/La)),
  16.        rotation=γ,cx=m,cy=n)
  17. }

実行結果は以下のようになります。

  1. r$>  u <- StandardEllipse(4,0,16,0,0,-64)
  2.      str(u)
  3. List of 5
  4.  $ a       : num 4
  5.  $ b       : num 2
  6.  $ rotation: num 0
  7.  $ cx      : num 0
  8.  $ cy      : num 0
  9. r$>  u <- StandardEllipse(16,0,4,0,0,-64)
  10.      str(u)
  11. List of 5
  12.  $ a       : num 4
  13.  $ b       : num 2
  14.  $ rotation: num 1.57
  15.  $ cx      : num 0
  16.  $ cy      : num 0
  17. r$>  u <- StandardEllipse(7,-6*sqrt(3),13,-2*sqrt(3),-34,-27)
  18.      str(u)
  19. List of 5
  20.  $ a       : num 4
  21.  $ b       : num 2
  22.  $ rotation: num 0.524
  23.  $ cx      : num 1.73
  24.  $ cy      : num 2
  25. r$>  u <- StandardEllipse(-7,6*sqrt(3),-13,2*sqrt(3),34,27)
  26.      str(u)
  27. List of 5
  28.  $ a       : num 4
  29.  $ b       : num 2
  30.  $ rotation: num 0.524
  31.  $ cx      : num 1.73
  32.  $ cy      : num 2
5.3. 縦長の標準形

作成した StandardEllipse関数は、式(4.11) によって a ≧ b となり、横長の標準形になります。x軸半径 a が短軸の場合は、横長の標準形を 90°傾けたと解釈します。

しかし、後述の行列の固有値を使った一般的な説明では、たいてい標準形を縦長に描いています。

まず、標準形の楕円を、

    (5.3.1)

と定めます。すると、a は、x軸半径です。普通のxy座標平面に描けば、a>b なら横長の楕円になります。b>a なら縦長の楕円です。

式(5.3.1) に  を乗じて、一般形にします。

                (5.3.2)

この式(5.3.2) は、傾いていないと考え、γ=0 です。係数は、 です。

式(4.1) によって描くと、係数A、B、C は回転(γ)によって下図のように変化します。

左の図は a>b のとき、右の図は b>a のときです。

横長の楕円の回転(a=4 ,  b=2)

縦長の楕円の回転(a=2 ,  b=4)

 上の左の図は、γ=0 のとき、係数Aは最小、係数Cは最大で、係数Bはゼロです。式は、

 

 上の右の図は、γ=0 のとき、Aが最大、Cが最小となり、式は、

 

です。

この2つの図は、左図のように、回転が 90°ズレた関係です。

横長の標準形を選択した、StandardEllipse関数は、

を、「X軸半径4、Y軸半径2の標準形の楕円を、反時計回りに π/6 回転し、X軸方向に√3、Y軸方向に2移動した楕円」と解釈します。

 縦長の標準形を採れば、2つの半径を入れ替えて、回転角度に90°を加算することになります。

この場合、与式は、「X軸半径2、Y軸半径4の標準形の楕円を、反時計回りに 2π/3 回転し、X軸方向に√3、Y軸方向に2移動した楕円」となります。

最初から、縦長の楕円を標準形とするなら、式(4.12) で回転していない式の係数A、C を計算します。

以下は変更がなく、

平方完成した式は、

その回転していない式は、

となります。

ので、 と求められます。

回転は、式(4.16) によることになって、

式(4.4) は、

となるので、

したがって、atan関数の性質から、 となります。

「3.楕円の諸元から一般形を求める」に挙げた maxima のスクリプトで30°刻みに式を求めると係数A、B、C は、以下のようになります。

γ

30°

60°

90°

120°

150°

180°

210°

240°

270°

γ’

30°

60°

90°

120°

150°

180°

A

4

7

13

16

13

7

4

7

13

16

B

0

-6√3

-6√3

0

6√3

6√3

0

-6√3

-6√3

0

C

16

13

7

4

7

13

16

13

7

4

u

12

12

12

12

12

12

12

12

12

12

-u-(C-A)

-24

-18

-6

0

-6

-18

-24

-18

-6

0

-u+(C-A)

0

-6

-18

-24

-18

-6

0

-6

-18

-24

tan γ’’

-Inf

-√3

-1/√3

0

1/√3

√3

-Inf

-√3

-1/√3

0

γ’’

-90°

-60°

-30°

0

30°

60°

-90°

-60°

-30°

0

γ’’+180°

90°

120°

150°

180°

ただし、B=0 のときは、A<C なら -Inf 、A>C なら 0 とします。

        

もし、 なら、

しかし、横長の標準形の場合と計算が異なるので、もう少し工夫が必要なようです。

5.4. 回転角の計算

回転角 γ の計算を考えて見ます

として、式(4.15) は、

式(4.16) は、

です。

u は、下図のように、楕円を回転しても変化しない定数となります。

式(4.15) と、式(4.16) のどちらかを決めるのは、γ=0 のときの係数A、C によります。γ=0 のときの標準形が横長(a>b)であれば、C>A となって、式(4.15) が適用されます。縦長(b>a)であれば、A>C となって、式(4.16) が適用されます。

式(4.16) は、

 

とすると、横長(a>b)であれば、

 

縦長(b>a)であれば、

 

と計算すれば、0から180°まで連続した回転角が表せます。

これで縦長の標準形を採った関数が作れます。StandardEllipse関数の必要個所を変更したVerticalEllipse関数を作りました。

  1. VerticalEllipse <- function(A,B,C,D,E,F){
  2.   if(A<0){A<- -A;B<- -B;C<- -C;D<- -D;E<- -E;F<- -F;}
  3.   if(abs(B)<24*.Machine$double.eps) {
  4.     if(C<A){ γ <- 0;     La <- C; Lb <- A;
  5.     } else { γ <- pi/2;  La <- A; Lb <- C; }
  6.   } else {
  7.     u <- sqrt(B^2+(C-A)^2); v<- C+A
  8.     La <- (-u+v)/2; Lb <- (u+v)/2
  9.     γ <- atan2(sqrt(u+(C-A)),sign(B)*sqrt(u-(C-A)))
  10.   }
  11.   g <- 4*A*C-B^2
  12.   m <- (B*E-2*C*D)/g
  13.   n <- (B*D-2*A*E)/g
  14.   Ua2b2 <- (A*E^2-B*D*E+C*D^2)/g - F
  15.   list(a=sqrt(abs(Ua2b2/Lb)),b=sqrt(abs(Ua2b2/La)),
  16.        rotation=γ,cx=m,cy=n)
  17. }

x軸半径 a=4、y軸半径 b=2 の楕円を、0°から30°刻みに270°まで回転して10個の楕円の式を作ります。それぞれ、StandardEllipse関数と比較してみます。

  1. g <- apply(rbind(seq(0,3*pi/2,by=pi/6)),1,
  2.            function(γ){GeneralFormsCoeffEllipse(4,2,γ,0,0)})[[1]]
  3. cat("----A----B--C--D--E---F  a--b-Ang--m--n  a--b-Ang--m--n\n")
  4. for(i in 1:length(g$A)){
  5.   u <- StandardEllipse(g$A[i],g$B[i],g$C[i],g$D[i],g$E[i],g$F[i])
  6.   v <- VerticalEllipse(g$A[i],g$B[i],g$C[i],g$D[i],g$E[i],g$F[i])
  7.   B <- round(g$B[i],7)
  8.   if(B==round(6*sqrt(3),7)) B <- " 6√3"
  9.   else if(B==round(-6*sqrt(3),7)) B <- "-6√3"
  10.   else B <- sprintf("%4.0f",B)
  11.   ur <- round(u$rotation/pi*180)
  12.   vr <- round(v$rotation/pi*180)
  13.   s1 <- sprintf("%2.0f %2.0f %s %2.0f %2.0f %2.0f %2.0f",
  14.                 i,g$A[i],B,g$C[i],g$D[i],g$E[i],g$F[i])
  15.   s2 <- sprintf("%2.0f %2.0f %3.0f %2.0f %2.0f",
  16.                 u$a,u$b,ur,u$cx,u$cy)
  17.   s3 <- sprintf("%2.0f %2.0f %3.0f %2.0f %2.0f",
  18.                 v$a,v$b,vr,v$cx,v$cy)
  19.   cat(s1,s2,s3,"\n")
  20. }

同じ式に対して、StandardEllipse関数はx軸半径 a=4、VerticalEllipse関数は a=2 となります。

4番目の式は、 で、StandardEllipse関数は90°回転と見なし、VerticalEllipse関数は、回転していないと判定します。

  1. ----A----B--C--D--E---F  a--b-Ang--m--n  a--b-Ang--m--n
  2.  1  4   -0 16  0  0 -64  4  2   0 -0 -0  2  4  90 -0 -0 
  3.  2  7 -63 13  0  0 -64  4  2  30 -0 -0  2  4 120 -0 -0 
  4.  3 13 -63  7  0  0 -64  4  2  60 -0 -0  2  4 150 -0 -0 
  5.  4 16   -0  4  0  0 -64  4  2  90 -0 -0  2  4   0 -0 -0 
  6.  5 13  63  7 -0 -0 -64  4  2 120  0  0  2  4  30  0  0 
  7.  6  7  63 13 -0 -0 -64  4  2 150  0  0  2  4  60  0  0 
  8.  7  4    0 16 -0 -0 -64  4  2   0  0  0  2  4  90  0  0 
  9.  8  7 -63 13  0  0 -64  4  2  30 -0 -0  2  4 120 -0 -0 
  10.  9 13 -63  7  0  0 -64  4  2  60 -0 -0  2  4 150 -0 -0 
  11. 10 16   -0  4  0  0 -64  4  2  90 -0 -0  2  4   0 -0 -0 
5.5. ゼロ判定と精度

VerticalEllipse関数では、係数Bのゼロ判定を以下のようにしました。

  1.   if(abs(B)<24*.Machine$double.eps) {
  2.   }

係数Bの計算が γ=0、90°、180°・・・に対して。 以下のような値を返すからです。

  1. r$> a <- 4; b <- 2; γ <- seq(0,3*pi/2,by=pi/2)
  2.     B <- 2*(b^2-a^2)*cos(γ)*sin(γ);
  3.     B
  4. [1]  0.000000e+00 -1.469528e-15  2.939055e-15 -4.408583e-15

90°の倍数の角度に対して、正弦、余弦は、一方が1、他方がゼロなので、係数Bはゼロとなります。

 なので、誤差は24倍になっています。

正弦と余弦の積の部分は、左図のようになります。

数値計算では、0° と 360° は同じではなく、0°から遠い方が誤差が大きくなっていきます。

不必要に大きな角度を計算しない方がよさそうです。

これは、0°以外は正確な値ではないからかも知れません。90°、360°は、π/2、2π を指定するわけで、定数πを演算して引数としています。

小さい値として、「.Machine$double.eps」を使いました。

  1. r$> .Machine$double.eps
  2. [1] 2.220446e-16
  3. r$> 2^(-52)
  4. [1] 2.220446e-16

64ビットの浮動小数点数の仮数部は52ビットのようです。1を表すと、仮数部の最上位のビットは 1/2 を表します。2 の -1乗です。最下位ビットは、 となると言うことのようです。

GRS80による地球の赤道半径は6,378,137mです。赤道の周の1mmは、中心角(経度)で、1.567856e-10 と計算されます。

これは、必要な精度の目安になります。

5.6. 定数項Fの役割

楕円、

                (5.6.1)

を例に採りましたが、定数項だけを変えた、

                (5.6.2)

を考えてみます。

前述のように、1つの楕円の式は、横長と縦長の2通りの標準形で説明できます。ここでは、横長の標準形だけを考えます。

式(4.11) から、

となります。

係数A、B、C から計算されるので、最初の例と変わりません。

  の計算も、係数Fに依存しないので変わりません。

変わるのは、 で、式(5.6.1) は、

式(5.6.2) は、

です。したがって、式(5.6.1) は、

となり、 となります。

式(5.6.2) は、

となり、 となります。

  1. r$> u <- StandardEllipse(7,-6*sqrt(3),13,-2*sqrt(3),-34,-12)
  2. r$> str(u)
  3. List of 5
  4.  $ a       : num 3.5
  5.  $ b       : num 1.75
  6.  $ rotation: num 0.524
  7.  $ cx      : num 1.73
  8.  $ cy      : num 2

式(5.3.1)と(5.3.2) は、いずれも    ですが、 は異なっています。

定数項Fの違いによって、それぞれ、 となりました。それぞれ、なので、U は、それぞれ、 です。

式(5.3.1) を、 の形の標準形にすれば、

        

式(5.3.2) は、

        

です。これを、

        

の形で表した標準形は、それぞれ、右辺の  で割って、

原点が楕円の中心で、平行移動していないときは、

        

です。定数項F は、式が、標準形を何倍したものかを示します。

平行移動した楕円の式の定数項は、平行移動によって生じた1次の項による定数を含んでいます。その値は回転の影響を受けています。

5.7. 平方完成

2つの式を考えます。

                (5.7.1)

                (5.7.2)

この2式の係数A、B、C、D、E は、同じです。定数項Fだけが異なっています。

係数Fに依存しない計算は共通です。

式(4.1) から、

を解いて、

式(4.5) から、

係数Fに依存する、 は、異なっていて、

式(5.7.1) は、

式(5.7.2) は、

です。

したがって、与式が傾いていない(γ=0)のときの式は、式(4.21) から、

与式(5.7.1) は、

                        (5.7.3)

展開すると、

                (5.7.4)

与式(5.7.2) は、

                        (5.7.5)

展開すると、

                (5.7.6)

です。

回転していない、平行移動した楕円の式は、以下のように平方完成できます。

        

 と置くと、        

        

定数項を右辺に移して   で割ると標準形になります。        

        

式(5.7.4)、(5.7.6) は、それぞれ与式に対して、平方完成する式です。平方完成する式は係数B=0で、楕円の軸と座標軸が平行です。

式(5.7.3)、(5.7.5) は、それぞれ与式に対して、平方完成した式です。

ここでは、最初に m、n を求めていますが、平方完成した式は、平行移動、あるいは楕円の中心座標を示します。

また、 と置いて、

式(5.7.3) から、

        

定数項を右辺に移して、定数項で割ると、

となり、x軸半径 a=4、y軸半径 b=2 と分かります。

式(5.7.5) から、

        

定数項を右辺に移して、定数項で割ると、

となり、x軸半径 a=7/2、y軸半径 b=7/4 と分かります。

5.8. 回転と係数の関係

平方完成するために、γ=0 のときの楕円の式を求めましたが、同様に、任意の γ について、係数を計算できます。

以下の、のケースを計算してみます。

楕円の形状は、楕円の2つの半径 a、b だけで決まります。座標平面に描く場合は、位置と向きが必要になります。

上の式は、すべて、a=4、b=2 の同じ楕円で、楕円の中心の座標は  です。向きだけが、異なっています。

係数A、C は、γ=0、90° で最小最大となる変化で、 の間の変化をします。係数B は、γ=0、90° でゼロとなります。B>0 は、γ が 90° 超であることを示します。

ただし、a>b の場合で、a<b の場合は、下段の図のように 90° ズレたグラフとなります。γ が90°超かどうかの判定も逆になります。

6. 標準形との点の対応
6.1. 標準形からグラフを描く

楕円の軸が、座標軸に重なっていない場合のグラフを描く方法として、標準形を求めるものと思います。

標準形の楕円を描いて、それを、回転、平行移動して、目的の楕円のグラフとします。

上の図は、赤丸の点における卯酉線を描くのに、xy平面に標準形の楕円を描き、回転、平行移動し、高さを求めて描いた様子です。楕円体の平面による断面の周は截線と言うようです。楕円体と平面の式から、zを消去すれば楕円の式になります。必要なのは、標準形を求める操作で求められた回転角度と、楕円の中心座標(cx、xy)だけです。まず、回転角度γの回転行列、

を適用して、cx、cy を、それぞれ、x、y に加算します。

まず、楕円体と、地点Aの緯度経度(φ、λ)から座標(Ax、Ay、Az)を求め赤い丸を描きます。

  1. b <- 0.6                # 楕円体の極半径
  2. DrawEarthEllipsoid(b=b,drawspheroid=FALSE,alpha=0.1,plane_aplha=0.8)
  3. φ <- pi/6; λ <- pi/4    # 地点Aの緯度、経度
  4. β <- atan(b*tan(φ))     # 更成緯度に
  5. Ax <- cos(β)*cos(λ); Ay <- cos(β)*sin(λ); Az <- b*sin(β)
  6. # 地点Aの座標の赤い丸
  7. spheres3d(Ax,Ay,Az,radius=0.02,col="red")
  8. # 更成緯度βの鉛直線とZ軸の交点
  9. q <- b*sin(β)-cos(β)*tan(φ)

q は、地点Aにおける鉛直方向の直線と地軸(z軸)との交点です。経度λの経線の作る平面と直交し、地点Aにおける鉛直方向の直線を含む平面が切断面です。楕円体を垂直な平面で切った垂直截線の一つを描きます。

楕円体の式と、平面の式から z を消去した2次曲線の式の係数を求めて、その係数から、楕円の標準形の諸元を求めます。

  1. R <- (Ay^2+Ax^2)/(q-Az)
  2. A <- Ax^2+b^2*R^2; B <- 2*Ax*Ay
  3. C <- Ay^2+b^2*R^2
  4. D <- -2*Ax*q*R; E <- -2*Ay*q*R
  5. F <- R^2*(q^2-b^2)
  6. ep <- StandardEllipse(A,B,C,D,E,F) # 標準形の諸元

楕円の諸元は以下のように算出されます。

  1. r$> str(ep)
  2. List of 5
  3.  $ a       : num 0.908
  4.  $ b       : num 0.654
  5.  $ rotation: num 2.36
  6.  $ cx      : num 0.206
  7.  $ cy      : num 0.206

この標準形の楕円の周上の点列は容易に算出できます。

  1. t <- seq(0,2*pi,length.out=361)
  2. x <- ep$a*cos(t); y <- ep$b*sin(t); z <- rep(0,length(x))
  3. lines3d(x,y,z,lwd=1,col="black")

これを、回転、平行移動すれば、楕円体の式と、平面の式から z を消去した2次曲線上の点列の座標になります。

  1. cr<-cos(ep$rotation); sr<-sin(ep$rotation)
  2. Q <- rbind(c(cr,-sr),c(sr,cr))
  3. g <- Q %*% rbind(x,y) # 回転
  4. x <- g[1,]+ep$cx; y <- g[2,]+ep$cy # 平行移動
  5. lines3d(x,y,z,lwd=1,col="orange")

赤道面に描かれた楕円の座標に相当する z を、切断面のへ面の式から計算すれば、一周分の垂直截線が描けます。

  1. z <- -(Ax*x+Ay*y-R*q)/R       # x,yからz算出
  2. lines3d(x,y,z,lwd=1,col="gray80") # 截線描画

しかし、垂直截線の一部を描くには、楕円体表面の座標から、標準形の座標の発生に使用する媒介変数の範囲を決める必要があります。

6.2. 一部を描く

全周を描く場合は対応は問題ではありませんが、上の図の赤丸Aから赤道との交点Pまでと言うように、一部を描くとします。すると、標準形のどこが赤丸A'で、どこが赤道との交点P'か知る必要が生じます。

先だって、切断面と赤道面の交線を引いて、交点を求めておきます。

  1. u <- Ax^2+Ay^2; v <- q-Az
  2. x <- c(-(Ay*sqrt(u*v^2-q^2*u^2)-Ax*q*u)/(u*v),
  3.        (Ay*sqrt(u*v^2-q^2*u^2)+Ax*q*u)/(u*v))
  4. y <- c((Ax*sqrt(u*v^2-q^2*u^2)+Ay*q*u)/(u*v),
  5.      -(Ax*sqrt(u*v^2-q^2*u^2)-Ay*q*u)/(u*v))
  6. spheres3d(x,y,0,radius=0.02,col="blue")
  7. DrawLine2p(c(x[1],y[1],0),c(x[2],y[2],0),lwd=1)

地点Aから赤道面に下した垂線と、赤道面との交点の座標は、(Ax、Ay、0)です。また、切断面と赤道の交点の一方は、(x[1]、y[1]) です。これを逆の変換をすれば良いので、cx、cy を、それぞれ、x、y から減じ、逆向きの回転行列を適用します。

  1. x <- c(Ax,x[1]) - ep$cx; y <- c(Ay,y[1]) - ep$cy
  2. g <- t(Q) %*% rbind(x,y) # 回転
  3. spheres3d(g[1,],g[2,],c(0,0),radius=0.02,col="red")

上の図の、A'、P' の赤丸が描かれます。これで、楕円体上の弧の描画範囲を、標準形の座標で取得できました。

変換行列は、以下です。

この楕円体上の弧の描画範囲に相当する、標準形の楕円の弧A'P'の座標があれば、全周を描いたのと同じ方法で、垂直截線の弧APを描くことができます。

6.3. 座標指定で標準形の弧を描く

標準形の楕円の周上の点を2つ指定して、その範囲の弧を描く点列を生成します。

楕円の周上を2点で区切れば、弧は2つできます。どちらが必要な方かを決める方法は思いつきません。

そこで、中心角の180°以下を短い方と考えて、長い方、短い方と呼び分けることにして、指定可能にします。作図結果を見て変えられるようにすることしかできません。しかも、ちょうど180°の場合は、どちらを指定すれば意図通りなのかは不明瞭です。

ここで、中心角は正の連続した値と考えています。しかし、atan()はーπ/2からπ/2、atan()は0からπ、ーπ/から0と言った値域の値を算出します。

そこで、座標から、中心角を0から2πで返す関数を作ります。x軸方向が0°です。

  1. a360 <- function(x,y, a,b){
  2.   u <- atan2(y, x);     u <- ifelse(u<0, 2*pi+u, u)
  3.   v <- atan2(a*y, b*x); v <- ifelse(v<0, 2*pi+v, v)
  4.   list(central=u,parametric=v)
  5. }

式(1,2) に従って、媒介変数となる角度も計算します。その上で、0°や360°を跨ぐ、連続した媒介変数を生成して、弧を表す点列の座標を求めます。

  1. SectionArc <- function(x1,y1, x2,y2, a,b,
  2.                         long=FALSE,
  3.                         col="black",lwd=1){
  4.   A1 <- a360(x1,y1, a,b); A2 <- a360(x2,y2, a,b)
  5.   nick <- pi/360
  6.   u <- A2$central - A1$central
  7.   if(u < 0) {u<- -u; tmp<-A1;A1<-A2;A2<-tmp}
  8.   a1 <- A1$parametric; a2 <- A2$parametric
  9.   if(((u<=pi) && (!long)) ||((u>pi) && (long))){
  10.     t <- seq(a1,a2,by=nick)
  11.   } else {
  12.     t <- c(seq(a2,2*pi,by=nick),
  13.            seq(0,a1,by=nick))
  14.   }
  15.   x <- a*cos(t); y <- b*sin(t)
  16.   invisible(list(x=x,y=y))
  17. }

この点列は、全周を描いたときと同じように描けます。まず、回転、平行移動します。

  1. cr<-cos(ep$rotation); sr<-sin(ep$rotation)
  2. Q <- rbind(c(cr,-sr),c(sr,cr))
  3. g <- Q %*% rbind(g$x,g$y) # 回転
  4. x <- g[1,]+ep$cx; y <- g[2,]+ep$cy
  5. lines3d(x,y,rep(0,length(x)),col="red",lwd=5)

zを計算して描けば、垂直截線の一部が描けます。

  1. z <- -(Ax*x+Ay*y-R*q)/R # x,yからz算出
  2. lines3d(x,y,z,lwd=2,col="blue") # 截線描画
6.4. 標準形から楕円体表面へ

上の図の、R、R、'R'' は、標準形のx軸方向の頂点Rの座標から、楕円体の表面の点R''を求めたものです。

これは、垂直截線の全周を描くのと同じです。

  1. g <- Q %*% rbind(u$a,0) # 回転
  2. x <- g[1,]+u$cx; y <- g[2,]+u$cy
  3. spheres3d(x,y,0,radius=0.02,col="black")
  4. z <- -(Ax*x+Ay*y-R*q)/R
  5. spheres3d(x,y,z,radius=0.02,col="black")

7. 行列の固有値との関係

2次曲線の一般形の式を、

と表し、係数行列Kの固有値、固有ベクトルを求めます。

与式が、

        (7.1)

のとき、

で、固有値、固有ベクトルは、与式の係数A、B、C だけから計算されます。

7.1. 固有値

与式(7.1) の係数行列 K の固有値、固有ベクトルは、Rのeigen関数で、

  1. r$> eigen(rbind(c(7,-3*sqrt(3)),c(-3*sqrt(3),13)))
  2. eigen() decomposition
  3. $values
  4. [1] 16  4
  5. $vectors
  6.            [,1]       [,2]
  7. [1,] -0.5000000 -0.8660254
  8. [2,]  0.8660254 -0.5000000

maximaのeigenvectors関数では、

固有値は、出力順にλ1、λ2とすると、いずれも、 です。

この値は、 式(5.1.1) と同じです。eigen関数がどのような計算をしているのかは分かりませんが、2x2の係数行列の固有値は、式(4.11) や式(4.12) で計算します。Rのeigen関数やmaximaのeigenvectors関数は、式(4.11) のように計算して、値の大きい方を先に出力するものと推測します。C++ のeigen3ライブラリでは、小さい値が先に出力されるので、出力順に決まりはないようです。

7.2. 固有ベクトル

Rのeigen関数の返した固有ベクトル P は、単位ベクトル化されています。

maximaの返した固有ベクトルは、単位ベクトル化できます。

固有ベクトルは2x2の行列で表されますが、列が2つのベクトルを表し、出力順で、2つの固有値に対応しています。

Rのeigen関数の固有ベクトルの第1列、固有ベクトル1は、 です。 です。

maximaの固有ベクトル1は、 となります。

また、単位ベクトル化したベクトルの要素は、三角関数と見なせます。

原点の周りに、反時計回りに γ 回転する回転行列は、

です。

Rのeigen関数の固有ベクトルは、120°回転する回転行列と同じです。

maxima の固有ベクトルは、-60°回転する回転行列です。

8. 固有値による標準形の計算

楕円の式、

                (8.1)

に対して、楕円の軸に平行な直線を座標軸としたx'y'座標を考えます。

x'y'座標の座標軸は、固有ベクトルによって与えられます。xy座標平面の原点を通るものをx'y'座標の軸とします。

x'y'座標で表した、与式の楕円の式、

                        (8.2)

を求めます。

その式は、固有値を  とすると、

        (8.3)

と表せます。

単位ベクトル化した固有ベクトル P は、

(8.4)

と、(x’、y’) から (x、y) への変換行列です。

また、

なので、

                                        (8.5)

です。

x'y'座標における平行移動、あるいは楕円の中心座標は、式(8.3) から

                        (8.6)        

です。

また、楕円の半径は、

                                (8.7)

です。

回転角 γ は、固有ベクトル P から求めます。

とすると、固有値 λ1 の固有ベクトルは、 です。

回転角 γ は、x'軸の傾きで、

        

です。

8.1. 係数A’、C’

固有値λ1、λ2 には、固有ベクトルPを介して以下の関係があります。

例えば、

の場合、

となります。

式(7.3.1) から、

となります。

楕円の軸に平行なx'軸、y'軸では、与式の係数A、B、C の項は、

                        (8.1.1)

です。

8.2. 係数D’、E’

与式の係数D、E の項は、式(8.4) から、

です。

したがって、

                        (8.2.1)

と変換できます。

9. eigen関数による計算例

楕円の式、

を例とします。

Rのeigen関数は、係数A、B、Cだけを使って以下のような結果を返します。

  1. r$> A<-7; B<- -6*sqrt(3); C<-13;
  2.     D<- -2*sqrt(3); E<- -34; F<- -27
  3. r$> M <- rbind(c(A,B/2),c(B/2,C))
  4. r$> eigen(M)
  5. eigen() decomposition
  6. $values
  7. [1] 16  4
  8. $vectors
  9.            [,1]       [,2]
  10. [1,] -0.5000000 -0.8660254
  11. [2,]  0.8660254 -0.5000000

eigen関数は、固有値として、16 と 4 、固有ベクトルとして、固有値16から  、固有値4から を出力します。

とインデックスが異なっています。

ベクトル は、互いに直交し、x'y'座標の座標軸を示します。xy座標の原点を通るものとします。

このベクトルは、楕円の軸に、それぞれ平行です。

式(8.5) から、

となるので、x'y'座標における楕円の式は、(8.1.1) と合わせて、

                        (9.1)

です。

式(8.6) は、        

となり、式(9.1) は、以下のように平方完成できます。

        

        

        

        

 とおくと、

        

定数を右辺に移して、定数で両辺を割ると、

        

        

したがって、与式が表す楕円は、x軸半径 a=2、Y軸半径 b=4 の楕円で、その中心は、x'y'座標の (m’、n’) です。

 (m’、n’) は、式(7.3.1) でxy座標に変換して、 となります。

回転角は、最初の固有ベクトル から、 となります。

10.  固有値による標準形の計算の図示

楕円の式、

は、x'y'座標平面で以下のように表されます。

                

これを、図示すると左図のようになります。

また、楕円の軸を座標軸X、Yとすれば、楕円は、

と表されます。

2つの固有ベクトルは互いに直交する直線を表し、X軸、Y軸に平行です。xy座標の原点を通る直線を x'軸、y'軸とします。

楕円の中心、あるいは、XY座標の原点、をx’y’座標で表して、 (m'、n') とします。

すると、上の2式は、それぞれ、x’y’座標で、  

と表されます。

平方完成した式、

は、

であることを示します。

与式の楕円は、x’y’座標で、  

と表され、この図を描いて、反時計回りに120°回転して、xy座標平面上の与式の楕円と重ねると下段の図となります。

11. 固有ベクトルによる座標変換

楕円の標準形を求めて、一般の楕円を以下のように描くことは、

① 媒介変数 θ によって単位円の周の座標値の列を作り、② 標準形の楕円にします。③ 固有ベクトル Pによって原点の周りに回転し、④ 平行移動すると解釈できます。

これは、すべてxy座標平面でのことです。

左図は、赤丸 (2、0) は固有ベクトルによって、青丸 (-1、√3) に変換されることを示しています。

この座標(2、0) 、(-1、√3) は、共にxy座標の座標値です。

これに対して、固有ベクトルが x'軸、y'軸を与えると考えると、図の青丸のx'y'座標で表した座標は、(2、0) で、そのxy座標での座標値は (-1、√3) だ、と言うこともできます。

これは、共に青丸の座標で、(2、0) は x'y'座標での座標値であり、(-1、√3) は xy座標での座標値です。

固有ベクトルは、xy座標で回転行列として機能します。

また、x'y'座標からxy座標へ、座標の値を変換します。

x'y'座標で表した楕円は、xy座標で表した回転していない楕円と同じです。

12. 行列による楕円の描画

楕円の媒介変数表示は、以下のように表せます。

XY座標は原点が楕円の中心で、楕円の軸が座標軸です。a はX軸半径、b はY軸半径です。

固有ベクトル P、楕円の中心座標 (m、n) が既知であるなら、任意の楕円は以下のようにして描画できます。

楕円の式が、

のとき、以下のスクリプトは、係数A、B、C、D、E、F から楕円を描きます。

  1. EllipseA <- function(A,B,C,D,E,F, col="black",lwd=1,lty=1){
  2.   Ua2b2 <- (A*E^2-B*D*E+C*D^2)/(4*A*C-B^2) - F
  3.   u <- eigen(rbind(c(A,B/2),c(B/2,C)))
  4.   λ1 <- u$values[1]; λ2 <- u$values[2]; P <- u$vectors
  5.   a <- sqrt(Ua2b2/λ1); b <- sqrt(Ua2b2/λ2)
  6.   DdEd <- (c(D,E) %*% P)[1,]
  7.   M <- P %*% rbind(c(a,0),c(0,b))
  8.   mn <- (P %*% c(-DdEd[1]/(2*λ1),-DdEd[2]/(2*λ2)))[,1]
  9.   t <- seq(0,2*pi,length.out=361)
  10.   g <- M %*% rbind(cos(t),sin(t)) + mn
  11.   lines(g[1,],g[2,],col=col,lwd=lwd)
  12. }

式(4.6) から、

係数A、B、C から求めた固有値をλ1、λ2、固有ベクトルを P とします。

 なので、

また、 なので、

とすると、

        

        

となり、

です。

XY座標において、座標 (-m'、-n') を原点とする XY座標に平行な x'y'座標を考えます。

XY座標の原点は、 x'y'座標で (m'、n') です。楕円の周上の点 (X、Y) は、 x'y'座標の座標値では (X+m'、Y+n') です。

また、 x'y'座標から xy座標への変換は、

なので、xy座標における楕円の中心は、

または、

        

        

と計算します。

13. 解法の逆順での楕円の描画

しかし、この手順は、標準形を求める手順の逆順ではありません。2x2の係数行列の固有値、固有ベクトルを使った解法の逆順は、

です。これは、

です。

楕円、

を、固有値、固有ベクトルを使った解法の逆順で描いてみます。

  1. plot(0,0,pch=20,xlim=c(-4,4),ylim=c(-6.5,1.5),
  2.      xlab="x'",ylab="y'")
  3. abline(h=0,v=0,col="blue")
  4. # 楕円の式の係数A,B,C,D,E,F
  5. A <- 7; B <- -6*sqrt(3); C <- 13
  6. D <- -2*sqrt(3); E <- -34; F <- -27
  7. # 係数A,B,Cから固有値、固有ベクトルを求める
  8. u <- eigen(rbind(c(A,B/2),c(B/2,C)))
  9. λ1 <- u$values[1]; λ2 <- u$values[2];
  10. P  <- u$vectors
  11. # 楕円の半径を計算
  12. Ua2b2 <- (A*E^2-B*D*E+C*D^2)/(4*A*C-B^2) - F
  13. a <- sqrt(Ua2b2/λ1); b <- sqrt(Ua2b2/λ2)
  14. # 係数D,Eからx'y'座標での係数D',E'を計算
  15. v <- (c(D,E) %*% P)
  16. Dd <- v[1,1]; Ed <- v[1,2]
  17. # x'y'座標での楕円の中心座標(m',n')
  18. md <- -Dd/(2*λ1); nd <- -Ed/(2*λ2)
  19. # X軸、Y軸をx'y'座標平面に描く
  20. abline(v=md,h=nd,col="cyan")
  21. # x'y'座標平面に楕円を描く
  22. t <- seq(0,2*pi,length.out=361)
  23. x <- a*cos(t)+md; y <- b*sin(t)+nd
  24. lines(x,y)

スクリプトの最後の行を、以下のようにすれば、下段の図のように、xy座標平面に描くことができます。

  1. g <- P %*% rbind(x,y)
  2. plot(g[1,],g[2,],t="l",xlim=c(-2,6),
  3.      ylim=c(-2,6),xlab="x",ylab="y")
  4. abline(v=0,h=0)
  5. h <- P %*% rbind(c(-8,8,0,0,md,md,-8,8),
  6.                  c(0,0,-8,8,-8,8,nd,nd))
  7. lines(h[1,1:2],h[2,1:2],col="blue")
  8. lines(h[1,3:4],h[2,3:4],col="blue")
  9. lines(h[1,5:6],h[2,5:6],col="cyan")
  10. lines(h[1,7:8],h[2,7:8],col="cyan")

楕円の半径は、固有値、固有ベクトルでは決まりません。しかし、 は、すべての係数A、B、C、D、E、F を使って計算できるので、固有値と合わせて、X軸半径 a。Y軸半径 b を求めることができます。

楕円は半径だけで決まるので、楕円の軸をX軸、Y軸として、直ちにXY座標平面に描くことができます。

14. 平方完成

楕円の式、

        

は、原点を中心とした楕円を描きます。これに対して、

        

は、中心が (m、n) の楕円を描きます。 平方完成は、 (m、n) を明らかにします。

これまでに、1つの与式について、2通りの平方完成した式、(14.1) と、(14.2) を挙げました。

式(4.1) は、楕円の式の係数が、半径 a、b と、 中心座標 m、n、および回転角 γ で決まることを示しています。

楕円の形状と位置は変わらず、回転角 γ だけを変えることを考えます。

楕円、

は、式(4.11) から、

        

 式(4.1) は、 γ =0 のとき、

となります。

与式は、γ =0 のとき、

        

と表せ、

        

と平方完成できます。

式(4.5) から  で、式(4.17) から です。

平方完成は、

                        (14.1)

で、平方完成をした元の式は、

です。

この場合、標準形はx軸半径 a=4 の横長で、反時計回りに30°回転したと解釈されます。

これに対して、固有値 λ1、λ2、固有ベクトル P による解法では、楕円の軸に平行な x'y'座標を考えて、

から、

        

を、平方完成します。平方完成した式は、

                        (14.2)

となり、 です。これは、x'y'座標で表した値です。

と、座標変換すると、 です。

15. 都合の良い固有値、固有ベクトル

eigen関数は、処理系によって異なる値を返します。固有値の出力順や、ベクトルの向きには決めごとがないためです。

Rのeigen関数で、以下の楕円の式を計算すると、x軸半径 2、y軸半径 4の楕円を120°回転して、x軸方向に√3、y軸方向に2 平行移動した平行移動した楕円と求められました。

この式は、「5. 楕円の諸元から一般形を求める」で、楕円の諸元を、

として作成したものです。

私にとって都合が良いのは、「x軸半径 4、y軸半径 2の楕円を30°回転して、・・・」となることです。

15.1. 固有値、固有ベクトルの計算

楕円の式の係数 A、B、C から固有値は式(4.11) によって以下のように計算できます。

係数Aが負なら式に-1を乗じて常に係数が正とすれば、

 

となります。式(4.12) によって、

        

と計算すれば、

 

となります。

固有ベクトルは、以下のように、 という値を採ることができます。

また、maxima の固有値は、以下のように解釈できます。

したがって、固有値 λ に対応する固有ベクトルは、(1、v) で、

          、または、

に対応した固有ベクトルは、楕円の軸に平行な x'軸の方向ベクトルです。x軸とx'軸の重なった状態が0°で、反時計回りにx'軸を回転します。180°で0°の楕円と重なります。こう考えると、固有ベクトルは、

が、都合が良く、

          、または、                        (15.1.1)

また、y'軸については、

から、

          、または、                (15.1.2)

15.2. 都合の良い計算

2x2の行列の固有値は2つの値で、その2つに特別な意味付けはありません。固有ベクトルも、いろいろに表現可能です

標準形の楕円は、x軸半径を a 、y軸半径を b として、

と表されます。分数を避けて、2次曲線の一般式に当てはめると、

係数行列の固有値は、

                 (7,11,1)

または、

                 (7,11,2)

となります。

式(4.11) で固有値を計算することに決めると、 と決まりますが、式(7,11,1) か、式(7,11,2) かは分かりません。

ここでは、式(7,11,2) と決めます。したがって、 となり、標準形を横長と考えることになります。

また、楕円の式は、

と表すことになります。定数項を右辺に移して、両辺を で割れば、

        

となり、 はx軸半径の2乗を表します。 に対応する固有ベクトルはx軸の方向ベクトルです。

楕円は180°回転すると、元の楕円と重なります。したがって、回転角 γ は、反時計回りに 0 ≦ γ < 180° で表すことにします。x軸と重なっていた楕円の軸が、反時計回りに回転し、x軸となす角が γ です。この楕円の軸を x’軸とすると、この軸のベクトルは1象限から2象限へ回転していきます。y'軸は2象限から3象限へ回転していきます。

したがって、 に対応する固有ベクトルは (x, 1) 、 に対応する固有ベクトルは (-1, y) ように求めれば良いことになります。

これをスクリプトにすると、

  1. StandardEllipseH <- function(A,B,C,D,E,F){
  2.   if(A<0){A<- -A;B<- -B;C<- -C;D<- -D;E<- -E;F<- -F;}
  3.   if(abs(B)<24*.Machine$double.eps) {
  4.     if(A<C){ λ1 <- A; λ2 <- C; P <- rbind(c(1,0),c(0,1))
  5.     } else { λ1 <- C; λ2 <- A; P <- rbind(c(0,-1),c(1,0)) }
  6.   } else {
  7.     u <- sqrt(B^2+(C-A)^2); v<- C+A
  8.     λ1 <- (-u+v)/2; λ2 <- (u+v)/2
  9.     v11 <- 2*(λ1-C)/B; v22 <- -2*(λ2-A)/B
  10.     r1 <- sqrt(v11^2+1); r2 <- sqrt(1+v22^2)
  11.     P <- rbind(c(v11/r1,-1/r2),
  12.                c(1/r1,  v22/r2))
  13.   }
  14.   γ <- atan2(P[2,1],P[1,1])
  15.   DdEd <- c(D,E) %*% P
  16.   m <- DdEd[1,1] / (2*λ1); n <- DdEd[1,2] / (2*λ2);
  17.   Ua2b2 = (λ1 * m^2) + (λ2 * n^2) - F
  18.   a = sqrt(Ua2b2 / λ1); b = sqrt(Ua2b2 / λ2);
  19.   v <- P %*% c(-m,-n)
  20.   list(a=a,b=b,angle=γ,cx=v[1,1],cy=v[2,1],vectors=P)
  21. }

まず、係数Aが負の場合は、式に-1を乗じて計算することにします。

係数Bがゼロの場合、楕円の軸は座標軸と平行です。この状態の回転が 0° なのか 90° なのかを係数A、Cで判断します。

回転していない横長の楕円は A<C になります。固有値は係数で、固有ベクトルは0°か90°かで決まり計算の必要はありません。

係数Bがゼロ以外の場合は、固有値、固有ベクトルを計算する必要があります。式(4.12) で固有値を求めます。

固有ベクトル P は、(x, 1) 、(-1, y) ように求めて、ベクトルのx軸に対する回転角が、それぞれ 1、2象限、2,3象限になるようにします。これを単位ベクトル化しておきます。

回転角 γ 、x軸半径 a、y軸半径 b を求めます。平行移動を表す m, n はx'y'座標で楕円の中心を表します。これを、xy座標に変換します。

  1. r$> u <- StandardEllipseH(
  2.     7,-6*sqrt(3),13,-2*sqrt(3),-34,-27)
  3.     str(u)
  4. List of 6
  5.  $ a      : num 4
  6.  $ b      : num 2
  7.  $ angle  : num 0.524
  8.  $ cx     : num 1.73
  9.  $ cy     : num 2
  10.  $ vectors: num [1:2, 1:2] 0.866 0.5 -0.5 0.866
  11. r$> u <- StandardEllipseH(4,0,16,0,0,-64)
  12.     str(u)
  13. List of 6
  14.  $ a      : num 4
  15.  $ b      : num 2
  16.  $ angle  : num 0
  17.  $ cx     : num 0
  18.  $ cy     : num 0
  19.  $ vectors: num [1:2, 1:2] 1 0 0 1
  20. r$> u <- StandardEllipseH(16,0,4,0,0,-64)
  21.     str(u)
  22. List of 6
  23.  $ a      : num 4
  24.  $ b      : num 2
  25.  $ angle  : num 1.57
  26.  $ cx     : num 0
  27.  $ cy     : num 0
  28.  $ vectors: num [1:2, 1:2] 0 1 -1 0

 

は、a=4、b=2 の楕円を反時計回りに30°回転した楕円をx軸方向に√3、y軸方向に2平行移動した楕円と求められます。

 

は、回転していないと判断され、

   

は、90°回転していると見なされます。

15.3. eigen関数の結果の組み替え

前述の StandardEllipseH() は、都合が良いように、固有値と固有ベクトルを計算しました。同様に、eigen関数の結果を組み替えることもできます。

係数行列 K に対して、

Rのeigen関数は、固有値、固有ベクトルを、

と出力し、縦長の標準形を120°回転したことを表しています。

横長の標準形になることを目的に、この結果の λ1 と λ2 を入れ替えます。

それに伴って、固有ベクトル P の列を入れ替えます。

楕円は180°の回転で元の楕円と重なることから、x軸と重なっていた楕円の軸が0から180°まで回転すると考えれば良いので、軸を表すベクトルは1象限、2象限の向きと考えます。y軸と重なっていた楕円の軸は、これに応じて、2象限、3象限を向きます。

固有ベクトル P の1列目の要素の符号を反転します。

これで、x'軸、y'軸の方向ベクトルは、それぞれ、 となります。

 は、30°です。また、

です。

 x'y'座標での与式、

は、

        

を、平方完成します。平方完成した式は、

        

です。x'y'座標での平行移動  です。xy座標では、

        

と変換して、 です。xxまた、

        

となります。

与式は、x軸半径 a=4、Y軸半径 b=2 の標準形の楕円を、反時計回りに30°回転し、x軸方向に√3 、y軸方向に 2 平行移動した楕円です。


題目一覧へmikeo_410@hotmail.com(updated: 2022/09/08)