2014-06-15

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

２変量の常微分方程式があって、２変数の時間変化が２つのベクトルを軸としてその２軸のそれぞれに指数関数の係数を与えた和で表されるとき、そのy=1平面への射影に複比保存が表れるのだが、それの「証明」というか、ひたすらな式変形で納得するためのメモ
常微分方程式
- ただし、この常微分方程式を表す $2\times 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(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$ によらないことである。それを式変形で確かめよう
まずの基本要素となるをと表すこととする
- $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}=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$ にもよらず、２つの固有値と $\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$
結局、２変量の常微分方程式の係数行列があったときその(実数)固有値の差 $\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$ となる
ちなみに、この複比を満足する点列は２つの固有値ベクトルの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