The order of accuracy for my Crank-Nicolson scheme is wrong.

Question

David Gillcrist 2021-7-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/884314-the-order-of-accuracy-for-my-crank-nicolson-scheme-is-wrong

I've been trying to figure out what the issue is in my numerical scheme I've written to solve the particle in the box problem for the Schrodinger equation. I've used a Crank-Nicolson scheme to implicitly solve for the wave function numerically, the issue I'm running into is that the Crank-Nicolson scheme is second order accurate in space, but I'm only finding that my numerical solution is on first order accurate. I am near 100% positive that the Crank-Nicolson scheme is order (2,2) for the Schrodinger equation, just like the heat equation, and that there should be no special exceptions here. Because of this, I'm assuming there's an error somewhere in my code for solving of the numerical solution.

My function for solving the particle in a box is here:

function [wave_function_final, wave_function, error] = Schrodinger_CrankNicolson_PIAB(h,k,T,xL,xR,N)
    m = ceil((xR - xL)/h); % Number of x-axis intervals
    h = (xR - xL)/m;
    n = ceil(T/k); % Number of t-axis intervals
    k = T/n;
    mu = k/h^2;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % natural units ---> MASS = 1, hbar = 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    DiagTerm = -2*(1 - 1i*2/mu);
    
    A = diag(DiagTerm*ones(1,m-1)) + diag(ones(1,m-2),1) + diag(ones(1,m-2),-1);
    
    B = -conj(A);
    
    wave_function = zeros(n+1,m-1);
    
    x_val = xL + h*(1:m-1);
    wave_function(1,:) = sqrt(2/(xR - xL))*sin(N*pi/(xR-xL) * (x_val - xL));
    
    for tstep = 2:(n+1)
        wave_function(tstep,:) = ((A'*A)\A'*B * wave_function(tstep-1,:)')';
    end
    
    wave_function = [zeros(n+1,1), wave_function, zeros(n+1,1)]; % <-- Concatenates Columns to the left and right of wave_function for boundary conditions
    
    wave_function_final = wave_function(end,:);
    
    En = (N*pi)^2/8;
    
    error = sqrt(sum( abs( exp(-1i*En*T)*wave_function(1,2:end-1) - wave_function_final(2:end-1) ).^2));
end

And my live script I use for running is here:

Ms = [5 10 20 40 80];
T = 16/pi;
xL = -1;
xR = 1;
N = 1;
orders = zeros(1,length(Ms)-3);
Errors = zeros(1,length(Ms));
xx = -1:0.05:1;
En = (N*pi)^2/8;
p = plot(xx,real(exp(-1i*En*T)*sqrt(2/(xR - xL))*sin(N*pi/(xR-xL) * (xx - xL)) ), 'Color', [0 0 0], 'LineStyle',"--", 'LineWidth',2);
p.Color(4) = 0.2;
for j = length(Ms):-1:1
    disp(j)
    h = 2/Ms(j);
    k = 1/Ms(j)^2;
    
    [wave_function_final, wave_function, error] = Schrodinger_CrankNicolson_PIAB(h,k,T,xL,xR,N);
    hold on
    xx = -1:h:1;
    plot(xx,real(wave_function_final),'LineWidth',1.5)
    Errors(j) = error;
end
for j = 4:length(Ms)
    orders(j-3) = log(abs( (Errors(j) - Errors(j-1))/(Errors(j-1) - Errors(j-2)) )) / log( abs( (Errors(j-1) - Errors(j-2))/(Errors(j-2) - Errors(j-3)) ) );
end
grid on
hold off
disp(Errors)
disp(orders)
figure 
plot(Errors)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

The order of accuracy for my Crank-Nicolson scheme is wrong.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

The order of accuracy for my Crank-Nicolson scheme is wrong.

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论