2変量常微分方程式と複比の関係

  • 2変量の常微分方程式があって、2変数の時間変化が2つのベクトルを軸としてその2軸のそれぞれに指数関数の係数を与えた和で表されるとき、そのy=1平面への射影に複比保存が表れるのだが、それの「証明」というか、ひたすらな式変形で納得するためのメモ
  • 常微分方程式\begin{pmatrix}\frac{dx}{dt}\\ \frac{dy}{dt} \end{pmatrix} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} k_1 & 0 \\ 0 & k_2 \end{pmatrix} \frac{1}{ad-bc} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix}
    • ただし、この常微分方程式を表す2\times 2行列は2つの実固有値を持つものとしてV \Sigma V^{-1}に分解してある(V=\begin{pmatrix} a & b \\ c & d \end{pmatrix},\Sigma =\begin{pmatrix} k_1 & 0 \\ 0 & k_2 \end{pmatrix},V^{-1} =  \frac{1}{ad-bc} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix})
    • この解は\begin{pmatrix} x(t) \\ y(t) \end{pmatrix} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} exp(k_1 \times t) & 0 \\ 0 & exp(k_2 \times t) \end{pmatrix} \frac{1}{ad-bc} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix} \begin{pmatrix} x_0 \\ y_0 \end{pmatrix}である
    • これを\begin{pmatrix} x(t) \\ y(t) \end{pmatrix} = V exp(\Sigma t) V^{-1} \begin{pmatrix} x_0 \\ y_0 \end{pmatrix}と書く
    • ここで、両辺に左から V^{-1}を掛けると
    • V^{-1}\begin{pmatrix} x(t) \\ y(t) \end{pmatrix} = V^{-1} V exp(\Sigma t) V^{-1} \begin{pmatrix} x_0 \\ y_0 \end{pmatrix}
    • V^{-1}\begin{pmatrix} x(t) \\ y(t) \end{pmatrix} =exp(\Sigma t) V^{-1} \begin{pmatrix} x_0 \\ y_0 \end{pmatrix}となるが、
  • ここで\begin{pmatrix} z \\ w \end{pmatrix} =V^{-1} \begin{pmatrix} x \\ y \end{pmatrix}と置き換える(\begin{pmatrix} x \\ y \end{pmatrix} = V \begin{pmatrix} z \\ w \end{pmatrix}と同じこと)
    • \begin{pmatrix} z(t) \\ w(t) \end{pmatrix} = exp(\Sigma t) V^{-1} \begin{pmatrix} z(0) \\ w(0) \end{pmatrix}となって、\Sigma^tが対角行列であるからz(t) = exp(k_1 t)z(0),w(t) = exp(k_2 t)w(0)と簡単に表せて
    • x(t) = a z(t) + b w(t),y(t)=cz(t) +d w(t)も単純に表せる
  • さて、t_1,t_2=t_1+\delta,t_3=t_1+2\delta,t_4=t_1+3\deltaという等間隔の4時刻を考える(x(t_i)x_iと簡略化して表すことにする)
  • y=1平面への原点から見た射影写像X_i=\frac{x_i}{y_i}とする
  • 今、このX_iの複比DR=\frac{(X_{3}-X_1)/(X_4-X_3)}{(X_2-X_1)/(X_4-X_2)}が一定である、というのはこの複比がt_1によらないことである。それを式変形で確かめよう
  • まずDRの基本要素となるX_u-X_vY_{u,v} = X_u-X_vと表すこととする
    • Y_{u,v}=\frac{x_u}{y_u}-\frac{x_v}{y_v}=\frac{x_uy_v-x_vy_u}{y_u y_v}である
    • これを使うとDR=\frac{Y_{3,1}/Y_{4,3}}{Y_{2,1}/Y_{4,2}}となり
    • さらに[tex:DR=\frac{*1/*2}{*3/*4}]
    • 気を付けると、y_iが消えるので
    • DR=\frac{(x_3y_1-x_1y_3)/(x_4y_3-x_3y_4)}{(x_2y_1-x_1y_2)/(x_4y_2-x_2y_4)}と簡単にできる
    • ここでさらにA_{i,j}=x_iy_j-x_jy_iと置き換えるとDR=\frac{A_{3,1}A_{4,2}}{A_{4,3}A_{2,1}}
  • A_{i,j}を指数関数で表そう
    • A_{i,j}=x_i y_j - x_j y_i
    • =(a z_0 exp(k_1t_i)+b w_0 exp(k_2t_i))(c z_0 exp(k_1t_j)+d w_0 exp(k_2t_j))
    •  - (a z_0 exp(k_1t_j)+b w_0 exp(k_2t_j))(c z_0 exp(k_1t_i)+d w_0 exp(k_2t_i))
  • よく注意するとac z_0^2 exp(k_1(t_i+t_j)),bd w_0^2 exp(k_2(t_i+t_j))の項は消え、
  • A_{i,j} = z_0w_0(ad-bc)(exp(k_1t_i+k_2t_j)-exp(k_1t_j+k_2t_i))となる
  • DR=\frac{A_{3,1}A_{4,2}}{A_{4,3}A_{2,1}}ではz_0w_0(ad-bc)は相殺できるから
  • DR=\frac{(exp(k_1t_3+k_2t_1)-exp(k_1t_1+k_2t_3))(exp(k_1t_4+k_2t_2)-exp(k_1t_2+k_2t_4))}{(exp(k_1t_4+k_2t_3)-exp(k_1t_3+k_2t_4))(exp(k_1t_2+k_2t_1)-exp(k_1t_1+k_2t_2))}
  • t_1,t_2=t_1+\delta,t_3=t_1+2\delta,t_4=t_1+3\deltaであるから、すべてのexp(k_1t_i+k_2t_j)t_i=t+(i-1)\deltaのうちのtの分は相殺されるから
  • DR=\frac{(exp(k_1(2\delta)+k_2(0\delta))-exp(k_1(0\delta)+k_2(2\delta)))(exp(k_1(3\delta)+k_2(\delta))-exp(k_1(\delta)+k_2(3\delta)))}{(exp(k_1(3\delta)+k_2(2\delta))-exp(k_1(2\delta)+k_2(3\delta)))(exp(k_1(\delta)+k_2(0\delta))-exp(k_1(0\delta)+k_2(\delta)))}
  • DR=\frac{(exp(2k_1\delta)-exp(2k_2\delta))(exp(3k_1\delta+k_2\delta)-exp(k_1\delta+3k_2\delta))}{(exp(3k_1\delta+2k_2\delta)-exp(2k_1\delta+3k_2\delta))(exp(k_1\delta)-exp(k_2\delta))}
  • 式はごちゃごちゃしているけれど、tは消えたので、どんなtをスタートにしても複比は一定である。また、z_0,w_0にもよらず、2つの固有値\deltaが複比の値を決めている
  • もう少し頑張ろう
  • [tex:DR=\frac{(exp(2k_1\delta)-exp(2k_2\delta))exp*5}{exp*6(exp(k_1\delta)-exp(k_2\delta))}]
  • DR=\frac{(exp(2k_1\delta)-exp(2k_2\delta))^2}{exp(k_1\delta+k_2\delta)(exp(k_1\delta)-exp(k_2\delta))^2}
  • exp(2k_1\delta)-exp(2k_2\delta)=(exp(k_1\delta)+exp(k_2\delta))(exp(k_1\delta)-exp(k_2\delta))に注意して
  • DR=\frac{(exp2k_1\delta)+exp(k_2\delta))^2}{exp(k_1\delta+k_2\delta)}
  • さらにDR=(exp(\frac{k_1-k_2}{2}\delta) + exp(\frac{k_2-k_1}{2}\delta))^2
  • 結局、2変量の常微分方程式の係数行列があったときその(実数)固有値の差\Delta_kと単位時間\deltaとによって、単位時間ごとの位置のy=1平面への原点を視点とする写像が作る点列は、複比\frac{\frac{X_3-X_1}{X_4-X_3}}{\frac{X_2-X_1}{X_4-X_2}}=(exp(\frac{\Delta_k}{2}\delta) + exp(\frac{-\Delta_k}{2}\delta))^2となる
  • ちなみに、この複比を満足する点列は2つの固有値ベクトルのy=1平面への写像を両端として、それらに収束する形になる
  • この式の素敵なところは、固有値の差になっていること。固有値複素数の場合も差にすれば、固有値は共役複素数なので複比は実数になってくる
  • やってみよう
n.iter <- 50
dr.calcs.min <- dr.calcs.max <- drs <- rep(0,n.iter)

for(ii in 1:n.iter){
	ab <- runif(2)
	ab[1]<-0.001*runif(1)
	lambdas <- c(ab[1]+1i*ab[2],ab[1]-1i*ab[2])
	aa <- rnorm(2)
	bb <- rnorm(2)
	vs <- cbind(aa+1i*bb,aa-1i*bb)
	A <- my.matrix.eigen(lambdas,vs)
	# 等間隔なt
	t <- seq(from=-20,to=20,length=100)
	X <- matrix(0,length(t),length(A[,1]))
	X.init <- runif(2)
	for(i in 1:length(t)){
	  X[i,] <- exp.m(A,t[i])[[1]] %*% X.init
	}
	x <- Re(X[,1]/X[,2])
	plot(Re(x))
	dr.calc <- my.doubleratio(x)

	delta <- t[2]-t[1]
	k1 <- lambdas[1]
	k2 <- lambdas[2]
	delta.k <- k1-k2
	dr <- (exp(delta.k/2*delta) + exp(-delta.k/2*delta))^2
	dr.calcs.min[ii] <- min(dr.calc)
	dr.calcs.max[ii] <- max(dr.calc)
	drs[ii] <- dr
}

plot(data.frame(dr.calcs.min,dr.calcs.max,drs))

*1:x_3y_1-x_1y_3)/(y_1y_3

*2:x_4y_3-x_3y_4)/(y_3y_4

*3:x_2y_1-x_1y_2)/(y_1y_2

*4:x_4y_2-x_2y_4)/(y_2y_4

*5:k_1+k_2)\delta)(exp(2k_1\delta)-exp(2k_2\delta

*6:2k_1+2k_2)\delta)(exp(k_1\delta)-exp(k_2\delta