熱方程式の数値解析(3)

なんてこった...。ダメもとでちょっと変えたらいけてしまった...。

本当はランダムウォークから熱方程式を導く場合には、
\begin{equation}
\frac{u(t + \Delta t, x) - u(t, x)}{\Delta t} = \frac{1}{2} \frac{u(t, x + \Delta x) + u(t, x - \Delta x) - 2 u(t, x)}{(\Delta x)^2}
\end{equation}

のように左右に遷移する確率が等しい=$\frac{1}{2}$に由来する係数が現れて、
\begin{equation}
\frac{(\Delta x)^2}{\Delta t} = D > 0
\end{equation}

の仮定で極限をとって、
\begin{equation}
\begin{cases}
u(0, x) = u_0(x), \\
\frac{\partial}{\partial t} u(t,\ x) = \frac{1}{2} D \Delta u(t,\ x), \quad t \ge 0,\ x \in \mathbb{R}^2
\end{cases}
\end{equation}

みたいになるようだが、係数なんてどうでもいいやろと適当にあしらったのがまずかったようだ...。
再び、$D=1$として、
\begin{equation}
u(t + \Delta t, x) - u(t, x) = \frac{1}{2}\big\{ u(t, x + \Delta x) + u(t, x - \Delta x) - 2 u(t, x) \big\}
\end{equation}

を差分方程式として数値解析させたらいい感じになった。

[initial_func.m]

% initial value at t=0

function y = initial_func(xx)
    length_xx = length(xx);
    y = zeros(1, length_xx);
    for i=1:length_xx
        x = xx(i);

        if x <= -1 || x >= 1
            y(i) = 0;
        else
            y(i) = exp(- 1 / (1 - x**2));
        end
    end
end

[Laplacian.m]

% Laplacian of uu

function y = Laplacian(uu)
    length_uu = length(uu);
    y = zeros(1, length_uu);

    for i = 1:length_uu
        u = uu(i);
        left_u  = 0;
        right_u = 0;
        if i - 1 >= 1
            left_u = uu(i - 1);
        end
        if i + 1 <= length_uu
            right_u = uu(i + 1);
        end

        y(i) = (left_u + right_u - 2*u) / 2;
    end
end

[heat.m]

function heat
    w  = 10; % half width of bar
    dx = 0.05;
    dt = dx**2;
    T  = 10; % total elapsed time [sec]

    n  = ceil(T/dt); % step num

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % sampling on x-axis
    xx = linspace(-w, w, ceil(2*w/dx));
    length_xx = length(xx);

    % solutions at each time
    u = zeros(n, length_xx);
    % initial value
    u(1,:) = initial_func(xx);

    % solve heat equation
    for step = 1:n-1
        t = dt * (step - 1); % current time
        u(step + 1,:) = u(step,:) + Laplacian(u(step,:));
    end

    show_result(u, xx, n);
    %show_result_animation(u, xx, n);
end

function show_result(u, xx, n)
    figure;
    hold on;
    plot(xx, u(1,:), "r");
    plot(xx, u(ceil(n/5),:), "b");
    plot(xx, u(ceil(2*n/5),:), "k");
    plot(xx, u(ceil(3*n/5),:), "g");
    plot(xx, u(ceil(4*n/5),:), "y");
    plot(xx, u(n,:), "m");
    axis([-10, 10, 0, 0.4]);
end

function show_result_animation(u, xx, n)
    figure;

    for i = 1:n
        plot(xx, u(i,:), "r");
        axis([-10, 10, 0, 0.4]);
        %sleep(0.01);
        refresh;
    end
end

show_resultで描かせた絵は以下。アニメーション版はこれが連続変形する。
熱がどんどん拡散していって、全体に温度が下がっていく様子が一応分かる。係数そんなに重要だったんか...。