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

常微分方程式の汎用解法（deSolve）

測地線を描く例にode関数を使ったものがありました。ODE(Ordinary Differential Equations)さえ知りませんでしたが、Rには、deSolve（General Solver for Ordinary Differential Equations）パッケージがあり、ode関数がありました。

測量の測地線は、GeographicLibで解説されているアルゴリズムで描いた測地線です。ODEによる線は、ode関数ごとRに写して描いてみました。少なくとも絵に描く分には完全に一致しています。

deSolveのode関数を使って同じように描ければ、測地線のアルゴリズムが理解できない私にとっては、ずっと見通しが良くなります。

1. deSolveのode関数の機能

微分方程式を解くとは、導関数から元になった未知関数を導くことのようです。

maximaで、2点確認しておきます。

１．傾き a が既知のとき、切片 b は、

\         b=y-ax

２．半径１の円の式 {x}^{2}+{y}^{2}=1 を微分すると、

\         2y\frac{dy}{dx}+2x=0

        したがって、

\         \frac{dy}{dx}=-\frac{x}{y}

導関数である微分方程式から、元の曲線を復元できるのは、各点の微分方程式の解である微分係数が接線の傾きを示すことによっているようです。

  1. 微分係数 <- function(x,y){-x/y}
  2. 切片 <- function(x,y){y-微分係数(x,y)*x}
  3. θ <- seq(pi/12,by=pi/6,length.out=12)
  4. x <- cos(θ); y <- sin(θ)
  5. for(i in 1:length(x))
  6. {
  7.   a <- 微分係数(x[i],y[i]); b <- 切片(x[i],y[i])
  8.   abline(b=a,a=b,col=rainbow(12)[i])
  9.   lines(c(0,x[i]),c(0,y[i]),col=rainbow(12)[i])
  10. }

導関数が \frac{dy}{dx}=-\frac{x}{y} であるとき、これを計算する関数 func(x,y) を以下のように書きます。

times は、－１から1まで、step刻みの数列で、算出される点の数と、点の間隔を決めます。ただし、y=0 は func() がゼロ除算となるので、両端は除きます。この例では、times は x です。

導関数から、半円が描かれました。

  1. func <- function(x,y,parms){list(c(-x/y))}
  2. step<-0.01
  3. times <- seq(-1+step,1-step,step)
  4. u <- ode(y=sin(acos(-1+step)), times, func, parms=NULL)
  5. plot(u[,1],u[,2],t="l")

ode関数の引数となる関数 func(x,y,parms) の parms は、未使用ですが記述が必要です。また、戻り値はリストである必要があり、リストには1つのベクトルが含まれます。

2. runge kutta法

測地線を、ode関数を利用して描く例がありましたが、Rのスクリプトではありません。引数もdeSolveパッケージのものとは異なっています。そこで、ode関数自体もRのスクリプトに写して最初の図を描きました。このode関数はデフォルトではrunge kutta法が選択されます。

  1. runge_kutta_ode <- function( f, y, range, step=0.001) {
  2.     xs <- seq(range[1], range[2]-step, by=step)
  3.     points <- matrix(ncol=length(xs)+1,nrow=length(y)+1)
  4.     points[,1] <- c(range[1],y);
  5.     row <- 2
  6.     for (x in xs) {
  7.       k1 <- f( x, y );
  8.       k2 <- f( x+step/2, y + k1*step/2 );
  9.       k3 <- f( x+step/2, y + k2*step/2 );
  10.       k4 <- f( x+step, y + k3*step );
  11.       y <- y + (k1 + 2*k2 + 2*k3 + k4)*step/6;
  12.       points[,row] <- c(x+step,y)
  13.       row <- row + 1
  14.     }
  15.     invisible(points);
  16. }

前述の半円を描いたdeSolveパッケージのode関数の使用例と同じことをしてみます。

  1. step <- 0.01; y0 <- sin(acos(-1+step));
  2. # deSolveのode
  3. funcA <- function(x,y,parms){list(c(-x/y))}
  4. times <- seq(-1+step,1-step,step)
  5. w <- ode(y=y0, times, funcA, parms=NULL)
  6. plot(w[,1],w[,2],t="l",col="cyan",lwd=5,
  7.      xlab="w[,1], u[1,]",ylab="w[,2], u[2,]")
  8. # ode_runge_kutta
  9. funcB <- function(x,y){-x/y}
  10. range <- c(-1+step,1-step)
  11. u <- runge_kutta_ode(f=funcB,y=y0,range=range,step)
  12. lines(u[1,],u[2,],col="navy")

2つのode関数は同じような半円を描きます。

  1. > range(u[1,]-w[,1])
  3. [1] -1.387778780781446e-17  2.220446049250313e-16
  4. > range(u[2,]-w[,2])
  5. [1] 0.000000000000000e+00 2.632459647861074e-05

しかし、一致するのは5桁程度のようです。

半径1の円ですから、(x,y)の2乗平均は1ですが、下段のグラフになりました。navyの色で描いた、runge_kutta_ode() の計算の方が、円に対する一致度が低いようです。

runge_kutta_ode() の計算を見てみます。step=0.1 として、x=-0.9 から -0.8 を求めます。

答えが分かっているので、初期値として、x=-0.9 に対応する、y0=sin(acos(-0.9)) を使います。

ユーザが定義する微分係数を計算する関数 f(x,y) を毎回4度呼び出します。この 4つの微分係数を k1、k2、k3、k4 とします。

        k1=f\left({{x}_{1}\ ,{y}_{1}}\right)=f\left({-0.9\ ,y0}\right)

        k2=f\left({{x}_{2}\ ,{y}_{2}}\right)=f\left({-0.9+\frac{step}{2}\ \ ,y0+k1\frac{step}{2}}\right)

        k3=f\left({{x}_{3}\ ,{y}_{3}}\right)=f\left({-0.9+\frac{step}{2}\ \ \ ,y0+k2\frac{step}{2}}\right)

        k4=f\left({{x}_{4}\ ,{y}_{4}}\right)=f\left({-0.9+step\ \ ,y0+k3\cdot{}step}\right)

関数 f(x,y) は、円の微分係数 -\frac{y}{x} を計算するので、

        k1=\frac{y0}{0.9}

        k2=\frac{y0+k1\left({step/2}\right)}{0.9+\left({step/2}\right)}

        :

点 \left({{x}_{1}\ ,{y}_{1}}\right)、\left({{x}_{2}\ ,{y}_{2}}\right)、\left({{x}_{3}\ ,{y}_{3}}\right)、\left({{x}_{4}\ ,{y}_{4}}\right) における微分係数 k1、k2、k3、k4 は、それぞれの点における接線の傾きで、図のようにxの区間で切った接線は、底辺の長さが同じ直角三角形と見なせます。

底辺の長さが同じなので、k1、2・k2、2・k3、k4 は各三角形の高さの比です。

これを平均したものを、y0に加算して、0.9+step=0.8 の y とします。

3. 複数の y

runge_kutta_ode() の計算は、引数 y が、Rのベクトルであっても同じように計算します。測地線を描く計算の y は 4つの値を持つベクトルになっています。ode関数から呼び出されるユーザが定義する関数は、ode関数の引数 y と同じ個数の値を返さなければなりません。

しかし、runge_kutta_ode() の処理の通り、yを構成する個々の値をode関数は見ていません。y を構成する値の意味はユーザが定義する関数だけの問題です。ode関数自体は、ユーザが定義する関数が返す値を複数の微分係数として、個々の y の値に適用するだけです。

この仕様がどんな応用を想定したものか分かりませんが、単純に並列に複数の y を計算してみます。

  1. rm(list=ls(all=TRUE))
  2. library(deSolve)
  3. source("lib/ode_runge_kutta.r")
  4. par(mar=c(3,3,2,1),mgp=c(1.5,0.5,0))
  5. # 2つの値を返す
  6. funcB <- function(x,y){c(-x/y[1],-x/y[2])}
  7. step <- 0.001; y0 <- c(1,-1);
  8. range <- c(0,1-step)
  9. u <- runge_kutta_ode(f=funcB,y=y0,range=range,step)
  10. plot(u[1,],u[2,],t="l",col="blue",ylim=c(-1,1),
  11.      ylab="u[3,],u[2,]")
  12. lines(u[1,],u[3,],col="orange")

0 から 0.9 までの範囲を指定して、yの初期値を1と-1にした2組のyを求めて半円を描きました。

引数 y に初期値2つを与えると、戻り値は3行になります。最初の行はｘ、2行目は1を初期値とした y、3行目は-1を初期値とした y になります。

4. ODEを使った測地線の描画の例

理解しようとしている例は、3軸の半径を指定する設定になっています。3軸方向の半径を、それぞれ a、b、c として、

        mean=\sqrt{\frac{{a}^{2}+{b}^{2}+{c}^{2}}{3}}

と計算しています。

ここでは赤道半径を1とした回転楕円体だけを考えます。

したがって、短軸半径、z軸半径、極半径 b は、扁平率を f として、

        b=1-f

です。また、離心率 e は、

        {e}^{2}=f\left({2-f}\right)

mean は、

        mean=\sqrt{\frac{2+{\left({1-f}\right)}^{2}}{3}}=\sqrt{1-\frac{{e}^{2}}{3}}

4.1. 座標

RGLで図を描くにはxyz直交座標の座標値を使います。

x=sinu\cdot{}cosv

y=sinu\cdot{}sinv

z=\left({1-f}\right)cosu

ただし、0\leq{}u\leq{}\pi{}\ ,\ \ 0\leq{}v\leq{}2\pi{}

例では楕円の表面を (u,v) で計算しており、 u は北極から南極へ向かって表した緯度、v は東周りに360°まで表した経度です。

uが示す角度は、楕円体上の点の中心角ではなく、更成緯度と同様に、左図のような補助球上の点を示します。

経度に相当する v は、赤道面であるxy平面の中心角ですが、cos(v) が x を指すので、x軸方向をゼロとして、東周りに2πまで測られます。

4.2. ode関数

deSolveパッケージのode関数にしたいのですが、前述の runge_kutta_ode() のようになっているので、この前提で調べます。

u0、v0 を楕円体上の位置、α0 を方位角の初期値としてode関数を呼び出しています。

  1. y0 <- c(u0, v0, cos(α0)*mean/sqrt(E(u0)), sin(α0)*mean/sin(u0))
  2. range <- c(0,2)
  3. u <- runge_kutta_ode(f=func, y=y0, range)

関数 E(u) は、

        E\left({u}\right)={cos}^{2}u+{\left({1-f}\right)}^{2}{sin}^{2}u=1-{e}^{2}{sin}^{2}u

4.3. ode関数の初期値

ode関数の初期値 y の4つの値を、\left({{u}_{0}\ ,\ {v}_{0}\ ,\ {g}_{0}\ ,\ {h}_{0}}\right)\ とすると、最初の2つは回転楕円体の表面の位置を表す緯度 u、経度 v の初期値です。

緯度 u は北極をゼロとしているので、常用の緯度を φ とすると

        \varphi{}=\frac{\pi{}}{2}-u 、u=\frac{\pi{}}{2}-\varphi{}

最初の図は、赤道上の経度ゼロから方位角45°で進むと言うものですが、{u}_{0}=\frac{\pi{}}{2}\ ,\ {v}_{0}=0 を初期値にしました。

方位角は、{\alpha{}}_{0}=\frac{\pi{}}{4}+\frac{\pi{}}{2} と、常用と90°異なった値で描いています。どう解釈するのか分かりませんが南極側から測れば、この角度になります。{g}_{0}\ ,\ {h}_{0} は、以下の値が初期値です。

        {g}_{0}=\frac{mean\cdot{}cos{\alpha{}}_{0}}{\sqrt{E\left({{u}_{0}}\right)}}=-1.04527154012371                (4.3.1)

        {h}_{0}=\frac{mean\cdot{}sin{\alpha{}}_{0}}{sin{u}_{0}}=0.6271629240742261                (4.3.2)

ode関数自体は y が1つの値でも、n個の値でも区別なく1つのベクトルとして同じ演算をします。戻り値は x を加えた 5行となります。x の行を付加する以外には、ode関数は値の意味を知りません。

ユーザが定義するode関数が呼び出す関数で意味を持ちます。

4.4. ode関数が呼び出す関数

ユーザが定義する、ode関数から呼び出される関数は、円を描く例では、x、y に対して微分係数を返す関数でした。

測地線を描く例のユーザ定義関数も同じくx、y の2つの引数を持ちます。ただし、使われているのは1つだけです。使われている引数 y は4つの値を持つRのベクトルです。

ユーザ定義関数の戻り値も4つの値を持つベクトルです。

この4つの組を \left({u\ ,v\ ,g\ ,h}\right)  と表すことにします。初期値は \left({{u}_{0}\ ,\ {v}_{0}\ ,\ {g}_{0}\ ,\ {h}_{0}}\right)\ で、ode関数の戻り値は、

        \begin{vmatrix} x_0 & x_1 & \cdot\cdot\cdot\\ u_0 & u_1 & \cdot\cdot\cdot\\ v_0 & v_1 & \cdot\cdot\cdot\\ g_0 & g_1 & \cdot\cdot\cdot\\ h_0 & h_1 & \cdot\cdot\cdot \end{vmatrix}

のような matrix です。

ode関数は、\left({{u}_{n}\ ,\ {v}_{n}\ ,\ {g}_{n}\ ,\ {h}_{n}}\right)\ から \left({{u}_{n+1}\ ,\ {v}_{n+1}\ ,\ {g}_{n+1}\ ,\ {h}_{n+1}}\right)\ を繰り返し計算するわけですが、ユーザ定義関数は、この間に複数回呼び出されます。

ユーザ定義関数は、引数 y の \left({u\ ,v\ ,g\ ,h}\right) から \left({u'\ ,v'\ ,g'\ ,h'}\right) を返します。

        u'=g

        v'=h

g'=\frac{\left({{f\left({2-f}\right)g}^{2}+{h}^{2}}\right)cosu\cdot{}sinu}{E\left({u}\right)}=\frac{\left({{e}^{2}{g}^{2}+{h}^{2}}\right)cosu\cdot{}sinu}{E\left({u}\right)}                 (4.4.1)

h'=\frac{-2gh}{tan\left({u}\right)}                                                 (4.4.2)

4.5. ode関数の処理

式(4.4.1)、(4.4.2) の計算を、

g\left({u,g,h}\right)=\frac{\left({{e}^{2}{g}^{2}+{h}^{2}}\right)cosu\cdot{}sinu}{E\left({u}\right)}                 (4.5.1)

h\left({u,g,h}\right)=\frac{-2gh}{tan\left({u}\right)}                         (4.5.2)

と記すことにします。

2点目をきめる runge_kutta_ode() の処理を記してみます。x に当たるユーザが定義関数の引数は省略して4つの値の組だけを書きます。

        k1=f\left({{u}_{0}\ ,\ {v}_{0}\ ,\ {g}_{0}\ ,\ {h}_{0}}\right)=\left({u'\ ,\ v'\ ,\ g'\ ,\ h'}\right)

        u'={g}_{0}

        v'={h}_{0}

        g'=g\left({{u}_{0}\ ,\ {g}_{0}\ ,\ {h}_{0}}\right)\

        h'=h\left({{u}_{0}\ ,\ {g}_{0}\ ,\ {h}_{0}}\right)

        k2=f\left({\left({{u}_{0}\ ,\ {v}_{0}\ ,\ {g}_{0}\ ,\ {h}_{0}}\right)+\left({u'\ ,\ v'\ ,\ g'\ ,\ h'}\right)\frac{step}{2}}\right)=\left({u''\ ,\ v''\ ,\ g''\ ,\ h''}\right)

        u''={u}_{0}\ +g'\frac{step}{2}

        v''={v}_{0}+h'\frac{step}{2}

        g''={g}_{0}+g\left({u'\ ,\ g'\ ,\ h'}\right)\ \frac{step}{2}

        h''={h}_{0}+h\left({u'\ ,\ g'\ ,\ h'}\right)\frac{step}{2}

        k3=f\left({\left({{u}_{0}\ ,\ {v}_{0}\ ,\ {g}_{0}\ ,\ {h}_{0}}\right)+\left({u''\ ,\ v''\ ,\ g''\ ,\ h''}\right)\frac{step}{2}}\right)=\left({u'''\ ,\ v'''\ ,\ g'''\ ,\ h'''}\right)

                ：

        k4=f\left({\left({{u}_{0}\ ,\ {v}_{0}\ ,\ {g}_{0}\ ,\ {h}_{0}}\right)+\left({u'''\ ,\ v'''\ ,\ g'''\ ,\ h'''}\right)}\right)=\left({u''''\ ,\ v''''\ ,\ g''''\ ,\ h''''}\right)

                ：

と計算して、

        {u}_{n+1}=\frac{u'+2u''+2u'''+u''''}{6}

        {v}_{n+1}=\frac{v'+2v''+2v'''+v''''}{6}

を求めることになります。これを maxima で行って、

        

のように長い式ができます。S は step です。長くてもmaximaからコピーしてRのスクリプトとして使えます。

しかし、見にくいので、g()、h() を先に計算して、

        {G_1}=g\left( {u_0},{g_0},{h_0}\right) \\ {H_1}=h\left( {u_0},{g_0},{h_0}\right)\\ {G_2}=g\left( \frac{2 {u_0}+S {g_0}}{2},\frac{2 {g_0}+{G_1} S}{2},\frac{2 {h_0}+{H_1} S}{2}\right)\\ {H_2}=h\left( \frac{2 {u_0}+S {g_0}}{2},\frac{2 {g_0}+{G_1} S}{2},\frac{2 {h_0}+{H_1} S}{2}\right)\\ {G_3}=g\left( \frac{4 {u_0}+2 S {g_0}+{G_1} {{S}^{2}}}{4},\frac{2 {g_0}+{G_2} S}{2},\frac{2 {h_0}+{H_2} S}{2}\right)\\ {H_3}=h\left( \frac{4 {u_0}+2 S {g_0}+{G_1} {{S}^{2}}}{4},\frac{2 {g_0}+{G_2} S}{2},\frac{2 {h_0}+{H_2} S}{2}\right)\\ {G_4}=g\left( \frac{2 {u_0}+2 S {g_0}+{G_2} {{S}^{2}}}{2},{g_0}+{G_3} S,{h_0}+{H_3} S\right)\\ {H_4}=h\left( \frac{2 {u_0}+2 S {g_0}+{G_2} {{S}^{2}}}{2},{g_0}+{G_3} S,{h_0}+{H_3} S\right)\\

１ステップ目は、

        \\ {u_{1}}=\frac{6 {u_0}+6 S {g_0}+{G_3} {{S}^{2}}+{G_2} {{S}^{2}}+{G_1} {{S}^{2}}}{6} \qquad(4.5.3)\\ {v_{1}}=\frac{6 {v_0}+6 S {h_0}+{H_3} {{S}^{2}}+{H_2} {{S}^{2}}+{H_1} {{S}^{2}}}{6} \qquad(4.5.4)\\ {g_{1}}=\frac{6 {g_0}+{G_4} S+2 {G_3} S+2 {G_2} S+{G_1} S}{6} \qquad(4.5.5)\\ {h_{1}}=\frac{6 {h_0}+{H_4} S+2 {H_3} S+2 {H_2} S+{H_1} S}{6} \qquad(4.5.6)\\

です。

この計算は、n=0 から １のステップを計算したものですが、以下のようにRのスクリプトを書いて測地線を描くことができます。

  1. f <- 0.4; b <- 1-f
  2. clear3d(); rgl.light();
  3. vectors3d(c(1.3,0,0),origin=c(-1.3,0,0),col="red"); text3d(1.4,0,0,"X",col="red")
  4. vectors3d(c(0,1.3,0),origin=c(0,-1.3,0),col="green"); text3d(0,1.4,0,"Y",col="green")
  5. vectors3d(c(0,0,0.78),origin=c(0,0,-0.78),col="blue"); text3d(0,0,0.88,"Z",col="blue")
  6. ee <- ellipsoid(c=b,col="white",alpha=0.4)  # 回転楕円体
  7. shade3d(ee)
  8. Circle3D(c(0,0,0), c(1,0,0), c(0,1,0),col="red",lwd=1) # 赤道
  9. sMax <- 2
  10. S <- 0.001
  11. N <- sMax / S + 1
  12. mean <- sqrt((2+(1-f)^2)/3)
  13. α0 <- pi/4+pi/2
  14. u0<-pi/2; v0<-0
  15. x <- function( u, v ) { sin(u) * cos(v) }
  16. y <- function( u, v ) { sin(u) * sin(v) }
  17. z <- function( u ) { (1-f) * cos(u) }
  18. E <- function( u ) { cos(u)^2 + (1-f)^2*sin(u)^2 }
  19. g <- function(u,g,h){ ((f*(2-f)*g^2+h^2)*cos(u)*sin(u))/E(u) }
  20. h <- function(u,g,h){ -(2*g*h)/tan(u) }
  21. g0 <- cos(α0)*mean/sqrt(E(u0)); h0 <- sin(α0)*mean/sin(u0)
  22. xs <- seq(0,sMax,by=0.001)
  23. u <- rep(NA,N); v <- rep(NA,N)
  24. gg <- rep(NA,N); hh <- rep(NA,N)
  25. u[1] <- u0; v[1] <- v0; gg[1] <- g0; hh[1] <- h0
  26. for(i in 1:(N-1)){
  27.   G1<-g(u[i],gg[i],hh[i])
  28.   H1<-h(u[i],gg[i],hh[i])
  29.   G2<-g((2*u[i]+S*gg[i])/2,(2*gg[i]+G1*S)/2,(2*hh[i]+H1*S)/2)
  30.   H2<-h((2*u[i]+S*gg[i])/2,(2*gg[i]+G1*S)/2,(2*hh[i]+H1*S)/2)
  31.   G3<-g((4*u[i]+2*S*gg[i]+G1*S^2)/4,(2*gg[i]+G2*S)/2,(2*hh[i]+H2*S)/2)
  32.   H3<-h((4*u[i]+2*S*gg[i]+G1*S^2)/4,(2*gg[i]+G2*S)/2,(2*hh[i]+H2*S)/2)
  33.   G4<-g((2*u[i]+2*S*gg[i]+G2*S^2)/2,gg[i]+G3*S,hh[i]+H3*S)
  34.   H4<-h((2*u[i]+2*S*gg[i]+G2*S^2)/2,gg[i]+G3*S,hh[i]+H3*S)
  35.   u[i+1] <- (6*u[i]+6*S*gg[i]+G3*S^2+G2*S^2+G1*S^2)/6
  36.   v[i+1] <- (6*v[i]+6*S*hh[i]+H3*S^2+H2*S^2+H1*S^2)/6
  37.   gg[i+1] <- (6*gg[i]+G4*S+2*G3*S+2*G2*S+G1*S)/6
  38.   hh[i+1] <- (6*hh[i]+H4*S+2*H3*S+2*H2*S+H1*S)/6
  39. }
  40. X <- x(u,v); Y <- y(u,v); Z <- z(u)
  41. lines3d(X,Y,Z,col="blue")
5. deSolveのode関数

RのdeSolveパッケージのode関数の引数は以下のようになっています。

ode(y, times, func, parms, method, …)

これは、測地線をode関数で描いている例の関数とは異なっています。x の範囲とステップで指定していたことが times の数列になっているようです。

引数 parms は、y と共にユーザが定義する関数に渡される値です。ユーザ定義関数は、必ず3つの引数を書くようになっていていますが、今回は未使用です。ユーザ定義関数は Rの list を返さなければ、その旨のエラーになります。list にはRのベクトルで y と同じ数の値を含めます。

測地線をode関数と、deSolveのode関数で測地線を計算し、u、v をそれぞれ比較してみます。

  1. # 共通
  2. f <- 0.4; b <- 1-f; e2 <- f*(2-f)
  3. mean <- sqrt( 1 - e2/3 )
  4. u0<-pi/2; v0<-0; α0 <- pi/4+pi/2
  5. vMax<-2; step<-0.001
  6. x <- function( u, v ) { sin(u) * cos(v) }
  7. y <- function( u, v ) { sin(u) * sin(v) }
  8. z <- function( u ) { (1-f)*cos(u) }
  9. E <- function( u ) { 1 - e2*sin(u)^2 }
  10. y0 <- c(u0, v0, cos(α0)*mean/sqrt(E(u0)), sin(α0)*mean/sin(u0))
  11. # runge_kutta_ode ---------------------------
  12. funcA <- function(x,y){
  13.   u <- y[1]; g <- y[3]; h <- y[4]
  14.   c( g, h,
  15.      ((f*(2-f)*g^2+h^2)*cos(u)*sin(u))/E(u),
  16.      -(2*g*h)/tan(u) )
  17. }
  18. w <- runge_kutta_ode(f=funcA,y=y0,range=c(0,vMax))
  19. # deSolve ode -------------------------------
  20. funcB <- function(time,y,parms){
  21.     list( funcA(time, y) )
  22. }
  23. times <- seq(0,vMax,by=0.001)
  24. q <- ode(y=y0,times=times,func=funcB,parms=NULL)

10桁程度の精度で一致します。

  1. > range(w[2,]-q[,2])
  2. [1] -7.930376e-10  4.596048e-10
  3. > range(w[3,]-q[,3])
  4. [1] -3.426163e-10  8.818457e-11
6. 知ったこと

常微分方程式は、「未知関数とその導関数からなる等式で定義される方程式である微分方程式の一種」のようです。

既知の導関数から未知関数を導くことが「解」のようです。

ode(x,y) は x に対応した y の数列を返し、曲線上の点列を解とします。

ode関数が呼び出すユーザが定義する関数は微分係数を返し、点(x,y)における曲線の接線の傾きを示すと考えられます。

測地線の計算のユーザ定義関数は、x に当たる引数を使っていません。step や times が示す増分は、ode関数が y に反映してユーザ定義関数に渡されることになります。

ユーザ定義関数は引数 y の4つの値 \left({u\ ,v\ ,g\ ,h}\right) だけを処理し、結果も同じ4つの値の組になります。

ode関数は汎用なので単に4つの独立した y を並列に計算しているだけだと考えられます。測地線の計算という目的はユーザ定義関数だけが知ることです。

測地線は4つの独立した曲線の点列を求める動作で求められ、4つのうち2つが、それぞれ緯度の変化 u と、経度の変化 v です。ただし、ユーザ定義関数は4つの値を関連付けて計算しており、微分係数を示すことになっている4つの戻り値は独立ではありません。

ユーザ定義関数は4つの値から微小区間における緯度、経度の変化を計算しています。この計算は回転楕円体の表面における慣性運動です。

ユーザ定義関数は引数の x に当たる値を使っていませんが、ode関数は x を増加させながら各ステップの計算をします。ユーザ定義関数が点(x,y)の微分係数を返すものとして、\left({{x}_{n}\ ,{y}_{n}}\right) の微分係数を求め、接線上を少し進んだ点\left({{x}_{n}+\Delta{}x\ ,{y}_{n}+\Delta{}y}\right)の微分係数を求め、\left({{x}_{n＋1}\ ,{y}_{n+1}}\right)に至ります。

このユーザ定義関数が使っていない x は、おそらく経度です。計算に必要なのは長さで、経度差に相当する長さは緯度で変わります。扁平率 f=0.4 の回転楕円体について、赤道上の経度ゼロから方位角45°で進む例を描きましたが、{u}_{0}=\frac{\pi{}}{2}\ ,\ {v}_{0}=0 と、式(4.3.1)、(4.3.2)が初期値でした。

初期値の計算の mean は、楕円体の3軸の軸長の2乗平均です。赤道半径が1の回転楕円体を考えます。球であれば mean は1で、扁平率が大きくなると小さくなっていきます。

g,h の初期値 g0、h0は、球であれば左図のようになります。u が北極から測られ、赤道がπ/2となるので、一般的な緯度の場合と正弦と余弦が逆になっています。

もし、始点が u0=π/2 でなければ、h0 は1/sin(u0)になります。

赤道上でsin(π/2)=1の長さに相当する経度の差を ｖ とします。緯度 u における同じ経度差 v の長さは sin(u0) になります。h0 を1/sin(u0) するのは緯度に依存しない長さを表すためです。

楕円体の場合、g0 は緯度 u に依存するE(u)で割られます。これは子午線弧長の計算の楕円積分の被積分関数の2乗です。

ode関数が x や times としているのは、ここでは進んだ長さだと考えられます。曲面を慣性運動する物体の進んだ距離を、step や timesの値の間隔ごとに計算していきます。ただし、この値は回転楕円体の扁平率に依存しています。

上の図は常用の方位角45°（π/4）で進んだときの、球面と、回転楕円体（扁平率0.4）の図です。

いすれも、times は 0 から π/2 です。

  1. > f <- 0.0; b <- 1-f; e2 <- f*(2-f)
  2. > vMax <- pi/2
  3. > times <- seq(0,vMax,length.out=round(vMax/0.001))
  4. > q <- ode(y=y0,times=times,func=funcC,parms=NULL)
  5. > N<-length(q[,2])
  6. > X<-x(q[N,2],q[N,3]);Y<-y(q[N,2],q[N,3]);Z<-z(q[N,2])
  7. > u <- XYZ2φλ(X,Y,Z,f)
  8. > X<-x(q[,2],q[,3]);Y<-y(q[,2],q[,3]);Z<-z(q[,2])
  9. > c(u$φ/pi*180,u$λ/pi*180)
  10. 45 90 
  11. > s = sum(sqrt((X[2:N]-X[1:(N-1)])^2+(Y[2:N]-Y[1:(N-1)])^2+(Z[2:N]-Z[1:(N-1)])^2))
  12. > c(s,s/pi*180)
  13. [1]  1.570796 89.999996

f=0 とした球の場合、到達点は常用の緯度、経度で、それぞれ 45°、90° です。xyz座標値から各点の距離を求めて積算したものも 90° です。半径が1の球なので中心角は弧の長さと同じ値です。

  1. > f <- 0.4; b <- 1-f; e2 <- f*(2-f)
  2. > vMax <- pi/2
  3. > times <- seq(0,vMax,length.out=round(vMax/0.001))
  4. > q <- ode(y=y0,times=times,func=funcC,parms=NULL)
  5. > N<-length(q[,2])
  6. > X<-x(q[N,2],q[N,3]);Y<-y(q[N,2],q[N,3]);Z<-z(q[N,2])
  7. > u <- XYZ2φλ(X,Y,Z,f)
  8. > X<-x(q[,2],q[,3]);Y<-y(q[,2],q[,3]);Z<-z(q[,2])
  9. > c(u$φ/pi*180,u$λ/pi*180)
  10. 42.67075 76.80985 
  11. > s = sum(sqrt((X[2:N]-X[1:(N-1)])^2+(Y[2:N]-Y[1:(N-1)])^2+(Z[2:N]-Z[1:(N-1)])^2))
  12. > c(s,s/pi*180)
  13. [1]  1.295312 74.215899

扁平率が f=0.4 の回転楕円体の場合、到達点は、およそ、42.67°、76.8° で球より小さい値になります。進んだ距離も74.2です。

7. 解について

測地線の計算は、direct() と inverse() と決まっていて、測地線の第1問題、第2問題などと説明されています。

測地線の計算は、測地線の曲線を表すようには説明されず、「方位角αでs進むと(φ、λ)に至る」か「(φ、λ)に至るには方位角αでs進む」となります。

2点間の距離は、inverse() で計算されます。2点間に測地線を引くには、inverse() によって、α、s を知り、[0,s] の間を適宜分割して direct() を繰り返して測地線を表す点列を得ていました。

ode関数を使った例は、測地線の曲線を表すもので、見事に測地線を表す点列が得られます。

しかし、direct() と inverse() のようには、まだ、理解が足りません。

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