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

まぁ、難しいもんだ。

[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)) / 10;
        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;
    end
end

[heat.m]

function heat
    w  = 2; % half width of bar
    dx = 0.05;
    dt = dx**2;
    T  = 1; % 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, dt);
    %show_result_animation(u, xx);
end

function show_result(u, xx, dt)
    figure;
    hold on;
    plot(xx, u(1,:), "r");
    plot(xx, u(ceil(0.01/dt),:), "b");
    plot(xx, u(ceil(0.02/dt),:), "g");
    axis([-2, 2, -0.02, 0.04]);
end

function show_result_animation(u, xx)
    figure;

    for i = 1:8
        plot(xx, u(i,:), "r");
        axis([-2, 2, -0.02, 0.04]);
        sleep(0.5);
        refresh;
    end
end

で実装してみたものの解がすぐにおかしくなってしまう。確かに熱方程式は非可逆の発展をするが、本来は平滑化効果で時間発展ですごく滑らかになるはずなのにかえってジグザグになってしまった。
実装がバグってるかもしれないが、そもそもこんな単純な実装ではダメなのかもしれない。ちょっと興味深い結果になってしまった。