追跡曲線

追跡曲線

　「追跡曲線」とインターネットで調べてみると、追跡曲線について様々な情報を得ることが可能である。サイトによっては微分方程式を解析的に解いていたり、図が掲載されているかもしれない。しかし、それは結果だけ示されているだけであり、「過程」が見えてこない。そこでここではある対象物に向かって進んでいく様子をコンピュータ、特にRstudioを使ってその軌跡がどのようなものを描くのかを解析に解くのではなく、数値的に解き、その動いている様子を動画にすることをここでの最終目標とする。

問題設定

　ここでは次のような問題を扱うことにしよう。 まず初めに、点\((1,0),(0,1),(-1,0),(0,-1) \ \)にそれぞれ \(A \ \),\(\ B \ \),\(\ C \ \),\(\ D \ \)がいて、 \(A \ \)は\(\ B \ \)へ、\(B \ \)は\(\ C \ \)へ、\(C \ \)は\(\ D \ \)へ、\(D \ \)は\(\ A \ \) へ向かってそれぞれ同じ速さ\(\ v \ \)で歩き出す。その時に描く曲線を求めることにしよう。 この曲線を解析的に解こうとすると連立微分方程式を解くことになる。そこでコンピュータの出番である。

追跡曲線を求める

　問題設定で与えられた曲線を数値的に求めるために離散的に考える。 分かりやすくするためにまず\(\ A \ \)が\(\ B \ \)へ向かって歩く現象だけを扱い、それを他の3つの現象にも適用することで曲線を描く方法をとる。

　時間を\(\ t \ \)とすると、\(A \ \)、\(B \ \)の座標は\(\ t\ \) の関数として描けるので、\(A \ \)の座標を\((x_1,y_1)\)、 \(B \ \)の座標を\((x_2,y_2)\)とすると、\(x_1=x_1(t),y_1=y_1(t),x_2=x_2(t),y_2=y_2(t)\) と書ける。ここで\(A \ \)の速度ベクトルを\(\ \mathbf{v} = (v_{x_1},v_{y_1}) \ \)とすると、ある時刻 \(t_1\)における速度ベクトル\(\ \textbf{v} \ \)は、丁寧に記述すると次のようになる。 \[\begin{eqnarray*} v_{x_1} &=& v・\frac{x_2(t_1)-x_1(t_1)}{\sqrt{(x_2(t_1)-x_1(t_1))^2+(y_2(t_1)-y_1(t_1))^2}} \\ v_{y_1} &=& v・\frac{y_2(t_1)-y_1(t_1)}{\sqrt{(x_2(t_1)-x_1(t_1))^2+(y_2(t_1)-y_1(t_1))^2}} \end{eqnarray*}\]

よって微分方程式を解くために4次のルンゲ・クッタ法を用いることにする。よって以下のコマンドを用意しておく。

RK4 <- function(dx.dt=function(t,x){}, h=1E-7, x0=1, t_start=0, t_end=1){  
                                   ##４次のルンゲ・クッタ法
  t <- seq(t_start, t_end, by=h)
  nsteps <- length(t)-1
  nvar <- length(x0)
  
  m <- (matrix(0,nrow=nsteps+1,ncol=nvar+1))
  m[1,] <- c(t_start,x0)
  m[,1] <- t
  for(i in 1:nsteps){
    x_n <- c(t(m[i,-1]))
    s_1 <- dx.dt(t[i],x_n)
    s_2 <- dx.dt(t[i]+h/2,x_n+h*s_1/2)
    s_3 <- dx.dt(t[i]+h/2,x_n+h*s_2/2)
    s_4 <- dx.dt(t[i]+h,x_n+h*s_3)
    m[i+1,-1] <- m[i,-1] + h/6*(s_1+2*s_2+2*s_3+s_4)
  }
  return(m)
}

このRK4を用いて問題設定での人数が4人の場合に限らず、一般のp人の場合にも適用できるようなpolygon関数を作成することを目標にする。 　上で示した数式をもとにそれぞれの速度ベクトル\(\ nxy\ \)を格納した関数\(\ ｇ\ \)を作り、そのあとは、初期値、

polygon <- function(p){                              ##polygon関数
  
  g <- function(t,x){ #曲線
  
    v=1    ##速度
  
    dx <- numeric(p)
    dy <- numeric(p)
    nxy <- numeric(2*p)
  
    for(i in 1:p){
      dx[i] <- x[2*i+1]-x[2*i-1]                     ##x方向に関する距離
      dy[i] <- x[2*(i+1)]-x[2*i]                     ##y方向に関する距離
      nxy[2*i-1] <- v*dx[i]/sqrt((dx[i]^2+dy[i]^2))  ##x方向の速度ベクトル
      nxy[2*i] <- v*dy[i]/sqrt((dx[i]^2+dy[i]^2))    ##y方向の速度ベクトル
    } 
    
    dx[p] <- x[1]-x[2*p-1]                           ##DとAとのx方向に関する距離
    dy[p] <- x[2]-x[2*p]                             ##DとAとのy方向に関する距離
    nxy[2*p-1] <- v*dx[p]/sqrt((dx[p]^2+dy[p]^2))    ##DとAのx方向の速度ベクトル
    nxy[2*p] <- v*dy[p]/sqrt((dx[p]^2+dy[p]^2))      ##DとAのy方向の速度ベクトル
  
    return(c(nxy))
  }

  ##初期値設定
                                  
  x <- function(p){                        #初期値の設定
    x <- numeric(2*p)
    for(i in 1:p){
      x[2*i-1] <- cos(2*(i-1)/p*pi)
      x[2*i] <- sin(2*(i-1)/p*pi )
    }
    return(x)
  }

  x0 <- c(x(p))                            #初期値
  

  h <- 10^(-4)                               #時間の刻み幅
  r <- RK4(g, h=h, x0, t_start=0, t_end=3)   #4次のルンゲ・クッタの引数の設定

  ##座標設定
  
  xlim <- c(-1,1)                            #xの範囲
  ylim <- c(-1,1)                            #yの範囲

  for(i in 1:p){                              ##グラフの重ね合わせ
  plot(r[,2*i],r[,2*i+1],xlab="",ylab="",type="l",xlim=xlim ,ylim=ylim)
  par(new=T)
  }

  for(i in 1:p){                              ##四角形
  segments(cos(2*(i-1)/p*pi),sin(2*(i-1)/p*pi),cos(2*i/p*pi),sin(2*i/p*pi))
  }
  g1 <<- g  
  r <<- r
  xlim <<- xlim
  ylim <<- ylim  
}

これでpolygon関数を定義することができた。このpolygon関数を用いることで、与えられた曲線を簡単に描くことができる。それを以下に図と動きの様子を表した動画の二つを併せて載せる。。

polygon(4)

len <- length(r[,1])
a <- floor(len/400)

for(j in 1:200){
  par(new=F)
  for(i in 1:4){
    plot(r[1:j*a,2*i],r[1:j*a,2*i+1],xlab="",ylab="",xlim=xlim,ylim=ylim,type="l")
    par(new=T)
  }
  p <- 4
  for(i in 1:p){                              ##四角形
  segments(cos(2*(i-1)/p*pi),sin(2*(i-1)/p*pi),cos(2*i/p*pi),sin(2*i/p*pi))
  }  
}

\(p=6,8\)の場合

　上のようにpolygon関数の

polygon(6)

len <- length(r[,1])
b <- floor(len/400)

for(j in 1:280){
  par(new=F)
  for(i in 1:6){
    plot(r[1:j*b,2*i],r[1:j*b,2*i+1],xlab="",ylab="",xlim=xlim,ylim=ylim,type="l")
    par(new=T)
  }
  p <- 6
  for(i in 1:p){                              ##六角形
  segments(cos(2*(i-1)/p*pi),sin(2*(i-1)/p*pi),cos(2*i/p*pi),sin(2*i/p*pi))
  }
}

polygon(8)

len <- length(r[,1])
c <- floor(len/400)
for(j in 1:350){
  par(new=F)
  for(i in 1:8){
    plot(r[1:j*c,2*i],r[1:j*c,2*i+1],xlab="",ylab="",xlim=xlim,ylim=ylim,type="l")
    par(new=T)
  }
  p <-8
  for(i in 1:p){                              ##八角形
  segments(cos(2*(i-1)/p*pi),sin(2*(i-1)/p*pi),cos(2*i/p*pi),sin(2*i/p*pi))
  }
}