タイムトンネルの数値計算

これを数値的に解いてみる。http://d.hatena.ne.jp/ita/20071218/p1
長さLの一次元系があり、両端が時間Tだけずれてつながっている。差分法で単純なオイラー法, f[x] +=dt*dfdt(f[])を考える。

拡散方程式

方程式は $\Large \frac{d}{dt}f(x,t)=\frac{\partial}{\partial x}^2 f(x,t)$ 。差分だと

#define XNUM 1000
#define dt 0.001
#define dx (L/XNUM)
double f[XNUM];

f[x] += dt* (2*f[x]-f[x-1]-f[x+1])/(dx*dx);

普通はt=0の初期値を設定して時間発展させる。しかし両端の時間がずれていると未来の値が必要になり計算ができない。そこでメッシュの大きさdx,dtを調節してx=0,t=0 とx=L, t=T を結ぶ線の上に初期値を設定し、そっから解を広げていくことを考える

□□□□□□　　　□□□□□□
□□□□□□　　　□□□□□■
□□□□□□　　　□□□□■□
□□□□□□　　　□□□■□□
□□□□□□　→　□□■□□□
□□□□□□　　　□■□□□□
■■■■■■　　　■□□□□□

通常は自分の値と左右の値から上の値を求めるが、

□○□
■■■

斜めになってるとこれができない。時間が一回微分、空間が二回微分になってるため。今、過去までの値が分かっていて、一段上を計算することを考える

□□□□■
□□Ａ■■
□Ｂ■■■
□■■■■
■■■■■

Ａのところを計算するには、一段下とその左右の値が必要。しかしＢはまだ決まっていない。また、どこか一箇所決めれば次々に連鎖反応で全部が決まるが、斜めに一周回ってきたときに矛盾が出ないためには初めの値は勝手に決められない。これらを一気に全部決めるには一段上の値すべてに関する連立一次方程式を立てて解けばいい。Ａ＝■−Ｂ−■という形で、■は既知の量。Ａ＋Ｂ＝■となる。正方行列で、対角要素と一つ上が１、それ以外０というもの。行列式０じゃなければ解ける。あとで考える。

追記

書いて布団入って一分で気がついた。解けない。隣接する２つを足すというのは移動平均みたいなもんで情報が失われるから逆変換はできない。各点の自由度に空間一回微分か時間一回微分の変数を追加して方程式を書き直すしかない。

空間軸をT/Lだけ傾けた座標系で書き直すと、fのx微分を $\Large f_x$ などと書くと
$\Large \frac{T^2}{L^2} f_{tt}=f_t +\frac{2T}{L}f_{tx} -f_{xx}$ 。速度 $\Large f_t=v$ とおけば

$\Large \frac{T^2}{L^2} v_t = v+\frac{2T}{L}v_x -f_{xx}$
$\Large f_t = v$

差分でやるなら

f[i] += dt * v[i];
v[i] += dt * (v[i] + 2*T/L* (v[i+1]-v[i-1])/(2*dx) - (2*f[i]-f[i+1]-f[i-1])/(dx*dx))*L*L/T/T;

シュレディンガー方程式バージョン

上とほぼ同じ。だけど流れを考えて確率が保存するような差分にしないとまずい。拡散は普通の差分で全濃度が保存するけど、シュレディンガー方程式は単純にやるとダメ。
波動関数を振幅と位相に分解し $\Large \phi(x) = r(x)\exp(i\theta(x))$ と置く。
$\Large i\frac{d}{dt}\phi=\frac{\partial}{\partial x}^2 \phi$ に代入すると、rのx微分を $\Large r_x$ などと書くと
$\Large i r_t -r \theta_t = 2 i r_x \theta_x + r \theta_{xx} + r_{xx} -r \theta_x^2$
実部と虚部を比較して $\Large r_t = 2r_x \theta_x + r \theta_{xx}$ および $\Large \theta_t = -r_{xx}/r + \theta_x^2$ 。
確率 $\Large R=r^2$ とおくと $\Large R_t = (r^2)_t =2 r r_t = 4r r_x \theta_x +2r^2 \theta_{xx} = (2r^2\theta_x)_x$ 。周期境界条件あるいは左右無限大でr=0なら $\Large \int_a^b dx R_t = [2r^2\theta_x]_a^b=0$ となりRの合計は保存。これを満たすような差分スキームが好ましい。まずrでなくRを使って方程式を書き直す。これを $\Large \theta_t = -r_{xx}/r +\theta_{xx}$ 、 $\Large R_x = 2 r r_x$ および $\Large R_{xx}=2r_x^2+2r r_{xx}$ を使って書き直すと、

$\Large \theta_t =R_{xx}/2R -R_x^2/4R^2 +\theta_x^2$ 。
$\Large R_t = (2R \theta_x)_x$

ここまで厳密。
離散化する場合、各格子点の中間で流れ $\Large 2R \theta_x$ を差分で $\Large J(i+1/2) = \frac{R(i)+R(i+1)}{2} \frac{\theta(i+1)-\theta(i)}{DX}$ と計算し、この分だけRが点iと点i+1の間で移動するとすれば全R合計は保存する。

int iL=(i+XN-1)%XN;
int iR=(i+1)%XN;
double thx = (th[i]-th[iL])/dx; // defined at i-1/2
double Rave = 0.5*(R[iL]+R[i]); // defined at i-1/2
J[i] = 2*Rave*thx; // defined at i-1/2
double Rx = (R[iR]-R[iL])/(2*dx); //defined at i
double Rxx = (2*R[i]-R[iR]-R[iL])/(dx*dx); //defined at i
double thxx = (2*th[i]-th[iR]-th[iL])/(dx*dx); //defined at i
R[i] += dt * (J[i]-J[i+1])/dx;
th[i] += dt * (Rxx/(2*R) - Rx*Rx/(4*R*R) +thxx);

拡散と同様に現在と左右の値から１ステップ後の値が出る。
任意次元三角形メッシュでやるなら、各ノード上で値が定義され、ノードを結ぶ辺の中点で流れが定義される。各ノードをボロノイ分割した多面体の面積と流れをかけて合計したものが増分となる。ここまでタイムトンネルのない場合。

ディラック方程式バージョン

$\Large \Psi_1 + i\frac{\partial}{\partial x}\Psi_4 = i \frac{\partial}{\partial t}\Psi_1$
$\Large -\Psi_4 + i\frac{\partial}{\partial x}\Psi_1 = i \frac{\partial}{\partial t}\Psi_4$
$\Large \Psi_2 + i\frac{\partial}{\partial x}\Psi_3 = i \frac{\partial}{\partial t}\Psi_2$
$\Large -\Psi_3 + i\frac{\partial}{\partial x}\Psi_2 = i \frac{\partial}{\partial t}\Psi_3$

を解く。差分だと

f1[x]+= dt*(f1[x]+I*(f4[x+1]-f4[x])/dx)/I;
f4[x]+= dt*(-f4[x]+I*(f1[x+1]-f1[x])/dx)/I;

とか。これなら順々に埋めていけるけど、左右対称でないのが気持ち悪い。いろいろ厳密に保存するには都合が悪そう。

前野さんのコメント考察

とりあえずタイムゲートが常設されてる場合は順次解いていけるようですね。そうじゃなくて一瞬だけゲートが開く場合は解けない感じです

普通の一次元周期系（両端の赤を接着）である時間だけゲートが開いて過去と未来がつながる場合、上の絵だと緑同士、紫同士を同一視するゲートがあると、これを実際つなげると右のようなトポロジーになります。ここで端から順次解いていくと筒が分岐する部分で解けなくなります。いままで別の点だったところがひっつくので、そこが同じ値になるように準備しておかないと駄目で、ある意味未来の情報が必要になるんで解けません。

というかこういう多様体の上で一般的に偏微分方程式を解く場合ってどうなるんだろう。