0
$\begingroup$

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?

$\endgroup$

marked as duplicate by Bill Greene, Community Jul 6 '17 at 12:25

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

Browse other questions tagged or ask your own question.