Numerical Solution of 1-D Heat Equation

확률론적인 방법으로 1-D heat equation을 푸는 방법을 간단하게 Mathematica로 구현해보았습니다. 아래의 그래프

는 $t = 0$ 일 때 $[-1, 1]$에서 $1$의 값을 갖고 나머지 지점에서는 모두 $0$인 초기조건이 주어졌을 때,

\begin{equation}\label{heat_eq} \frac{\partial u}{\partial t} = \frac{1}{2}\Delta u \end{equation}

를 Monte-Carlo method로 계산하여 $(t, x) \in [0, 1] \times [-2, 2]$ 영역에서 밀도그래프로 그린 것입니다. (왜 굳이 \eqref{heat_eq}에서 $\Delta$ 앞에 $\frac{1}{2}$를 붙였냐면, Brownian motion의 infinitesimal generator가 $\frac{1}{2}\Delta$이기 때문입니다.) 이를 계산하기 위하여 5,000개의 Brownian sample path

를 생성하여 이용하였습니다. 참고로 코드는 다음과 같습니다.

(* Description.
u  Initial condition. Compiled for better performance.
xrng  Spatial variable x. {From, To, # of subdivision}.
trng  Time variable t. {From, To, # of subdivision}.
numSampB  # of Brownian sample paths.
*)
u = Compile[{{x, _Real}}, UnitStep[x + 1] - UnitStep[x - 1]];
xrng = {-2, 2, 50};
trng = {0, 1, 30};
numSampB = 5000;

Timing[
   data = {};
   dx = (xrng[[2]] - xrng[[1]])/xrng[[3]];
   dt = (trng[[2]] - trng[[1]])/trng[[3]];
   sqdt = Sqrt[dt];
   sampX = Table[xrng[[1]] + dx*i, {i, 0, xrng[[3]]}];
   sampB = ConstantArray[0, numSampB];
   For[ti = 0, trng[[3]] >= ti,
    If[ti > 0, 
     sampB += RandomVariate[NormalDistribution[0, sqdt], numSampB];];
    dataT = {};
    For[xi = 1, Length[sampX] >= xi, 
     AppendTo[dataT, Mean[u /@ (sampX[[xi]] - sampB)]]; xi++;];
    AppendTo[data, dataT];
    ti++;
    ];
   Clear[dx, dt, sqdt, sampX, sampB, dataT, ti, xi];
   ][[1]];
Print[StringForm["Time Elapsed: ``(s)", Round[%, 10^-4] // N]];
ListDensityPlot[data, ColorFunction -> "Temperature", PlotRange -> All]
ListPlot3D[data, ColorFunction -> "Temperature", PlotRange -> All]

확률론적 방법이 기존의 방법에 대하여 갖는 장점은 역시 초기조건의 regularity에 의존해서 알고리즘의 속도가 현저히 바뀌거나 하지 않는다는 점이지요. 단점이라면, 무식하게 Monte-Carlo method로 푸는 경우 샘플이 적으면 그래프가 안 예쁩니다.