# Comparison of different methods for solving $\partial_tA=\nabla^2_xA$ [duplicate]

For testing I would like to solve the equation $$\partial_tA=\nabla^2_xA,\,A_0=\sin(x)$$ In order to solve it, I use an implicit solver $$\left(1-\frac{\nabla^2_x}{2}dt\right)A_{n+1} = \left(1+\frac{\nabla^2_x}{2}dt\right)A_n$$ or the Runge-Kutta-method in two different orders. For solving the whole equation I use the following matlab script:

clear all;
warning('off', 'all');
x_max = 1;
x_num = 1000;
a0_val = 1;
b0_val = 1;
dx = x_max/20000;
steps = 1/dx;

x_vec = linspace(0, x_max, x_num).';
A0_vec = sin(x_vec);%linspace(1, x_num, x_num).';
ones_vec = ones(x_num, 1);
dh = (abs(x_vec(2)-x_vec(1)));
nabla_matrix = spdiags([-1*ones_vec 16*ones_vec -30*ones_vec 16*ones_vec -1*ones_vec], -2:2, x_num, x_num);
nabla_matrix(1, 1) = 35;
nabla_matrix(1, 2) = -104;
nabla_matrix(1, 3) = 114;
nabla_matrix(1, 4) = -56;
nabla_matrix(1, 5) = 11;
nabla_matrix(2, 1) = 11;
nabla_matrix(2, 2) = -20;
nabla_matrix(2, 3) = 6;
nabla_matrix(2, 4) = 4;
nabla_matrix(2, 5) = -1;
nabla_matrix(x_num, x_num) = 35;
nabla_matrix(x_num, x_num-1) = -104;
nabla_matrix(x_num, x_num-2) = 114;
nabla_matrix(x_num, x_num-3) = -56;
nabla_matrix(x_num, x_num-4) = 11;
nabla_matrix(x_num-1, x_num) = 11;
nabla_matrix(x_num-1, x_num-1) = -20;
nabla_matrix(x_num-1, x_num-2) = 6;
nabla_matrix(x_num-1, x_num-3) = 4;
nabla_matrix(x_num-1, x_num-4) = -1;
nabla_matrix = nabla_matrix / (12*(dh*dh));

A_equation = @(A_vec) ((nabla_matrix*A_vec.'));

A_matrix_single = speye(x_num)-dx/2.*nabla_matrix;
B_matrix_single = speye(x_num)+dx/2.*nabla_matrix;

A_matrix = zeros(steps, x_num);
A_matrix_NL_rk4 = zeros(steps, x_num);
A_matrix_euler = zeros(steps, x_num);
A_matrix_rk4 = zeros(steps, x_num);
A_matrix_rk42 = zeros(steps, x_num);

A_matrix(1,:) = A0_vec;
A_matrix_NL_rk4(1,:) = A0_vec;
A_matrix_euler(1,:) = A0_vec;
A_matrix_rk4(1,:) = A0_vec;
A_matrix_rk42(1,:) = A0_vec;

%Do steps
for i = 2:steps
A_k1 = A_equation(A_matrix_rk4(i-1,:));
A_k2 = A_equation(A_matrix_rk4(i-1,:)+A_k1.'*dx/2);
A_k3 = A_equation(A_matrix_rk4(i-1,:)-A_k1.'*dx+2*dx*A_k2.');
A_matrix_rk4(i,:) = A_matrix_rk4(i-1,:)+dx*(A_k1/6+2/3*A_k2+A_k3/6).';
end

for i = 2:steps
A_k1 = A_equation(A_matrix_rk42(i-1,:));
A_k2 = A_equation(A_matrix_rk42(i-1,:)+A_k1.'*dx/2);
A_k3 = A_equation(A_matrix_rk42(i-1,:)+0.5*dx*A_k2.');
A_k4 = A_equation(A_matrix_rk42(i-1,:)+dx*A_k3.');
A_matrix_rk42(i,:) = A_matrix_rk42(i-1,:)+dx*(A_k1/6+1/3*A_k2+A_k3/3+A_k4/6).';
end

for i = 2:steps
A_k1 = A_equation(A_matrix_NL_rk4(i-1,:));
A_k2 = A_equation(A_matrix_NL_rk4(i-1,:)+A_k1.'*dx/2);
A_k3 = A_equation(A_matrix_NL_rk4(i-1,:)+0.5*dx*A_k2.');
A_k4 = A_equation(A_matrix_NL_rk4(i-1,:)+dx*A_k3.');
A_matrix_NL_rk4(i,:) = speye(x_num)\(speye(x_num)*A_matrix_NL_rk4(i-1,:).'+speye(x_num)*(dx*(A_k1/6+1/3*A_k2+A_k3/3+A_k4/6)));
end

for i = 2:steps
A_matrix(i,:) = A_matrix_single\(B_matrix_single*A_matrix(i-1,:).');
end

figure;
imagesc(A_matrix);
title('Forward step')
figure;
imagesc(A_matrix_rk4);
title('RK4 step')
figure;
imagesc(A_matrix_rk42);
title('RK42 step')
figure;
imagesc(A_matrix_NL_rk4);
title('RK4_NL step')


but the results are way of my expectations. I get a more or less correct result when using the implicit approach, but the latter three all explode.

Why is that? Is my general approach flawed?

• This question has been asked many times here in one form or another, e.g. scicomp.stackexchange.com/questions/549/…. Beyond that, there are many references that are easily available by doing a simple search, e.g. scholarpedia.org/article/Method_of_lines. Jul 6 '17 at 12:17
• @BillGreene: I did not see that this question is equivalent to the heat equation, thus solving my problem. Thanks! Jul 6 '17 at 12:25