解を出したい

  • こちらで2次元座標に周期的状態の相表示がされている
  • ロトカ=ヴォルテラのごくごく簡単な係数表現はx+y-\log(x)-\log(y)=Cである
  • この式の固定点は(1,1)である
  • 一方、ぐるりと回るものの典型に「円」がある
  • 原点を中心とした円とロトカ=ヴォルテラの周回曲線との関係を考えたい
  • (x-\log(x)-1)+(y-\log(y)-1)=R^2(0 < x,y \le \infty)とz^2+w^2=R^2(-\infty \le x,w \le \infty)とが対応すると決める
  • 前式ではf(x)=x-\log(x)-1、後式ではf(z)=z^2と考えれば、どちらもf(x)+f(y)=R^2という曲線であることがわかる
  • ということは、x-\log(x)z^2との間をつなぐ関数があれば正円とロトカ=ヴォルテラの周回曲線とが対応付けできる
  • 今、正円はz,wについて対称、ロトカ=ヴォルテラも今のそれはx,yについて対称なので、x=y上の点に関する対応関係を考える
  • 2(x-\log(x)-1)=R^2,2(z^2)=R^2なので、2(x-\log(x)-1)=R^2の2つの解と2z^2=R^2の2つの解を相互に対応づけることとしよう
  • さて、そこから先が今回の記事の本題
  • 2(x-\log(x)-1)=R^2の解がほしいが、解けないので、計算機にお願いすることにする
  • RのrootSolveパッケージを使う。Newton-Raphsonで求める
    • 単一の等式 fun=0の解を求めるのがuniroot.all()関数
library(rootSolve)
R<-0.8
fun <- function (x) x-log(x)-1-R^2
All <- uniroot.all(fun,c(0,100),maxiter=10000)
curve(fun(x),0,2*max(All),main="uniroot.all")
points(All,y=rep(0,length(All)),pch=16,cex=2)
    • 連立方程式の解を求めるのがmultiroot()関数
      • 以下の例は、x^2+y^2=1とx=2yの交点座標を出している(
model <- function(x) c(F1= x[1]^2+x[2]^2-1,
F2= x[1]-2*x[2])
# first solution
(ss<-multiroot(model,c(1,1),useFortran=FALSE))
(ss<-multiroot(f=model,start=c(1,0)))
# second solution; use different start values
(ss<-multiroot(model,c(-1,0)))