mikeo_410
  1. ３．R
    1. １．Rの常套句
    2. ２．いろいろ（R）
    3. ３．座標の引数（ベクトル、matrix、bind、list）
    4. ４．plot
    5. ５．RGL
    6. ６．ラグランジュの未定乗数法(Rsolnp)
    7. ７．常微分方程式の汎用解法（deSolve）

ラグランジュの未定乗数法(Rsolnp)

円の周上の点P(x,y) から水平、垂直に直線を引いて楕円に内接する長方形を描くことを考えます。

円の半径を1として考えます。

点Pは長方形の頂点で、P(1,0) のとき面積はゼロです。

P\left({\frac{1}{\sqrt{2}}\ ,\ \frac{1}{\sqrt{2}}}\right) のとき、面積は 2 です。面積Sは、

        S\left({x,y}\right)=\left|{2x}\right|\cdot{}\left|{2y}\right|=4\left|{xy}\right|          (1)

と計算できます。

ただし、x、y は、円周上の点に限られ、以下の式を満たします。        

        G\left({x,y}\right)={x}^{2}+{y}^{2}=1           (2)

円周上の点は媒介変数 θ を使って生成することができます。

        x=cos\theta{}\ ,\ \ y=sin\theta{}

図の点Pが第1象限にあるときのxと面積の関係は左図のように描けます。

  1. > θ <- seq(0,pi/2,length.out=91)
  2. > x <- cos(θ); y <- sin(θ)
  3. > S <-  function(x,y){4*abs(x*y)}
  4. > par(mar=c(3,3,1.5,1),mgp=c(1.5, 0.5, 0))
  5. > plot(x,S(x,y),t="l",col="navy",cex.main=0.9,
  6. +      main="円に内接する長方形の面積(第1象限)")
  7. > (i <- which.max(S(x,y)))
  8. [1] 46
  9. > S(x[i],y[i])
  10. [1] 2
  11. > abline(v=x[i],h=S(x[i],y[i]),col="gray60")

媒介変数 θ が 1°刻みになっています。Rの配列のインデクスは1から数えるので、46番目は45°に相当します。この時、長方形の面積は2です。範囲を絞って、刻みを細かくすれば、面積の最大値を数桁程度計算することは十分できます。

ラグランジュの未定乗数法は、未定乗数 λ を使って、拘束条件がある関数の極値を求める方法のようです。

Rには、Rsolnpパッケージに solnp() があります。

式(1) が主関数、式(2)を等式制約関数と考えます。

  1. > library(Rsolnp)
  2. > 主関数 <- function(X){x<-X[1]; y<-X[2]; 4*abs(x*y)}
  3. > 等式制約関数 <- function(X){x<-X[1]; y<-X[2]; x^2+y^2}
  4. > u <- solnp(pars = c(1,1), fun=主関数,
  5. +            eqfun=等式制約関数, eqB=1)
  7. Iter: 1 fn: 2.2500       Pars:  0.75000 0.75000
  8. Iter: 2 fn: 2.0069       Pars:  0.70833 0.70833
  9. Iter: 3 fn: 2.0000       Pars:  0.70711 0.70711
  10. Iter: 4 fn: 2.0000       Pars:  0.70711 0.70711
  11. Iter: 5 fn: 2.0000       Pars:  0.70711 0.70711
  12. solnp--> Completed in 5 iterations
  13. > str(u)
  14. List of 10
  15.  $ pars       : num [1:2] 0.707 0.707
  16.  $ convergence: num 0
  17.  $ values     : num [1:6] 4 2.25 2.01 2 2 ...
  18.  $ lagrange   : num [1, 1] 2
  19.  $ hessian    : num [1:2, 1:2] 1 0 0 1
  20.  $ ineqx0     : NULL
  21.  $ nfuneval   : num 50
  22.  $ outer.iter : num 5
  23.  $ elapsed    : 'difftime' num 0.0508050918579102
  24.   ..- attr(*, "units")= chr "secs"
  25.  $ vscale     : num [1:4] 2.00 1.00e-08 7.07e-01 7.07e-01
  26. > c(u$pars[1],u$pars[2],1/sqrt(2))
  27. [1] 0.7071067843267312 0.7071067780463615 0.7071067811865475
  28. > u$pars-1/sqrt(2)
  29. [1]  3.140183779670735e-09 -3.140186000116785e-09

eqB=1 は、等式制約関数が返す値です。等式制約関数は、1つのベクトルを返し、複数の値を持つことができるようになっているようです。等式制約が複数ある場合、等式制約の式の数だけ値をベクトルで返し、それぞれ eqB にも同じ数の等式制約の式の右辺に当たる値を書きます。

また、主関数、等式制約関数の引数は1つのベクトルで、複数の値を含みます。引数がいくつあるかは pars によって認識されるものと思います。

pars は、主関数の引数の数だけ値を指定するようですが、どのように決めるのか良く分かりません。

結果は8桁程度の精度のようです。

精度指定は可能のようです。

  1. > ε <- .Machine$double.eps*64
  2. > u <- solnp(pars = c(1,1), fun=主関数,
  3. +            eqfun=等式制約関数, eqB=1,
  4. +            LB=c(ε,ε))
  6. Iter: 1 fn: 2.2500       Pars:  0.75000 0.75000
  7. Iter: 2 fn: 2.0069       Pars:  0.70833 0.70833
  8. Iter: 3 fn: 2.0000       Pars:  0.70711 0.70711
  9. Iter: 4 fn: 2.0000       Pars:  0.70711 0.70711
  10. Iter: 5 fn: 2.0000       Pars:  0.70711 0.70711
  11. solnp--> Completed in 5 iterations
  12. > u$pars-1/sqrt(2)
  13. [1] 1.110223024625157e-16 1.110223024625157e-16

この計算は、以下のように、連立方程式を作って、解くことのようです。

主関数 f\left({x,y}\right)=4xy 、等式制約関数 g\left({x,y}\right)={x}^{2}+{y}^{2} を定めます。

等式制約関数に未定乗数 λ を乗じて主関数から引いた式を作ります。

        F\left({x,y}\right)=f\left({x,y}\right)-\lambda{}\cdot{}g\left({x,y}\right)=4xy-\lambda{}\left({{x}^{2}+{y}^{2}}\right)

この式を、xで微分、yで微分して、2式を作り、等式制約関数と合わせた3式を連立方程式として解きます。

測地線の説明にある方位角

測地線の説明を私は理解できません。まず、α0、α1、α2 の方位角が何を指すか明確にしたいと思います。

赤道上の1点から、方位角 α0=α1=45° で出発して、出発点からの経度差 60° の地点の緯度 φ2 を求めます。

測地線の計算は、2点間の測地線長を求める invers() と、方位角αでs進むとどこに達するかを求める direct() の2つの関数で出来ています。

まず、direct() で、経度差 λ=60° となる距離sを求めます。

  1. > a <- 6378137
  2. > f <- 1/298.257222101; b <- 1-f
  3. > φ1 <- 0; β1 <- φ2β(φ1)
  4. > φ2 <- pi/6; β2 <- φ2β(φ2)
  5. > s <- seq(0,pi,length.out=181)
  6. > x <- rep(NA,length(s)); y <- rep(NA,length(s)); z <- rep(NA,length(s))
  7. > φ <- rep(NA,length(s)); λ <- rep(NA,length(s)); α <- rep(NA,length(s))
  8. > β <- rep(NA,length(s))
  9. > for(i in 1:length(s)){
  10. +     u <- GE_DirectA(0,0,pi/4,s[i],f)
  11. +     φ[i] <- u$φ2; λ[i] <- u$λ2; α[i] <- u$α2
  12. +     β[i] <- u$β2
  13. +     v <- φλ2XYZ(u$φ2, u$λ2, f)
  14. +     x[i] <- v$x; y[i] <- v$y; z[i] <- v$z
  15. + }
  16. > λ[67:70]/pi*180
  17. [1] 57.87028487208667 59.09296428544720 60.32899947471504 61.57820473966889  
  18. > λ[69]/pi*180
  19. [1] 60.32899947471504
  20. > φ[69]/pi*180
  21. [1] 41.12754870570981
  22. > s[69]
  23. [1] 1.186823891356144
  24. > s[69]*a
  25. [1] 7569725.373942603

等間隔で発生した距離 s の69番目が経度差 60° に最も近く、そのときの s[69] の値を使います。

  1. > main_func <- function(X) {
  2. +   s <- X[1];
  3. +   u <- GE_DirectA(0,0,pi/4,s,f)
  4. +   abs(pi/3-u$λ2)
  5. + }
  6. > eqfun <-  function(X) {
  7. +   s <- X[1];
  8. +   u1 <- GE_DirectA(0,0,pi/4,s,f)
  9. +   u2 <- GE_InverseA(0,0,u1$φ2,u1$λ2, f)
  10. +   u2$s - s
  11. + }
  12. > ε <- .Machine$double.eps*64
  13. > s <- solnp(s[69],main_func,eqfun=eqfun,LB=ε)$pars
  14. Iter: 1 fn: 0.000000001506       Pars:  1.18220
  15. Iter: 2 fn: 0.000000001506       Pars:  1.18220
  16. solnp--> Completed in 2 iterations
  17. > s*a
  18. [1] 7540211.394253859
  19. > u <- GE_DirectA(0,0,pi/4,s,f)
  20. > u$φ2/pi*180
  21. [1] 41.03451556619601
  22. > u$β2/pi*180
  23. [1] 40.93924575376187
  24. > u$λ2/pi*180
  25. [1] 59.99999991373227
  26. > u$α2/pi*180
  27. [1] 69.40062520997895
  28. > u$α21/pi*180
  29. [1] 249.400625209979

s[69] の値を使って、経度差が60°となる距離sを探します。

求めた s を使って、φ2、λ2を求めます。

この値を、国土地理院の測量計算サイトに入力します。

測地線長は、7,540,211.394(m)、方位角は、出発点→到着点が45°、到着点→出発点が249.400625°と表示されます。

プログラムと数式の関数が同じに見えない私にはピンときませんが、プログラムの関数が返す値が引数に応じて凸や凹であればなんでも極値を知ることができるようです。

しかし、測地線の説明自体に Lagrangian function L が使われているので、理解している人にとってはものすごく冗長なことをしているのだろうと思います。

測量計算サイトが距離としている測地線長の計算を理解したい訳ですが、赤道から方位角45°で経度差60°まで進んだときの到達点に関して、手元の計算結果と矛盾はないようです。

赤道上の点における法線は地心に向かうので方位角は補助球でも同じ角度になります。経度差60°まで進んだ点では、回転楕円体と、その補助球では異なる方位角が算出されます。

  1. > f <- 1/298.257222101; b <- 1-f
  2. > φ2 <- 41.03451556619601 /180*pi # 到達点の緯度、経度
  3. > λ2 <- 59.99999991373227 /180*pi
  4. > β2 <- φ2β(φ2)                   # 到達点の更成緯度
  5. > # 始点のxyz座標
  6. > p0 <- φλ2XYZ(0,0, f)
  7. > # 経度λ2の経線の作る平面
  8. > P60 <- Plane(c(0,0,0),c(0,0,1),c(cos(λ2),sin(λ2),0))
  9. > # 回転楕円体での方位角α21
  10. > p2e <- φλ2XYZ(φ2,λ2 , f)        # 到達点のxyz座標
  11. > z <- p2e$z*((b^2-1)/b^2)        # 法線と自転軸の交点
  12. > Pe <- Plane(c(0,0,z),p0,p2e)
  13. > 360 - Angle2Planes(P60,Pe)/pi*180
  14. [1] 249.4336724562589
  15. > # 補助球表面での方位角α21
  16. > p2s <- φλ2XYZ(β2,λ2 , 0)
  17. > Ps <- Plane(c(0,0,0),p0,p2s)
  18. > 360 - Angle2Planes(P60,Ps)/pi*180
  19. [1] 249.2413740967018

しかし、いずれも測量計算サイトの値にはなりません。まだ理解が足りないようです。

題目一覧へmikeo_410@hotmail.com(updated: 2023/01/21)