Error when Using Crank-Nicolson method to solve Schrodinger Equation

41 views (last 30 days)
Hello! I would like to know where I did wrong in my code, I am a fresh user of MATLAB and I just dipped into coding recently, MATLAB did not show exact error:
function [x, t, psi, psire, psiim, psimod, prob, v] = ...
sch_1d_cn(tmax, level, lambda, idtype, idpar, vtype, vpar)
% Inputs
%
% tmax: Maximum integration time
% level: Discretization level
% lambda: dt/dx
% idtype: Selects initial data type
% idpar: Vector of initial data parameters
% vtype: Selects potential type
% vpar: Vector of potential parameters
%
% Outputs
%
% x: Vector of x coordinates [nx]
% t: Vector of t coordinates [nt]
% psi: Array of computed psi values [nt x nx]
% psire Array of computed psi_re values [nt x nx]
% psiim Array of computed psi_im values [nt x nx]
% psimod Array of computed sqrt(psi psi*) values [nt x nx]
% prob Array of computed running integral values [nt x nx]
% v Array of potential values [nx]
% Define mesh and derived parameters ...
nx = 2^level + 1;
x = linspace(0.0, 1.0, nx)
dx = 2^(-(level));
dt = lambda * dx;
nt = round(tmax / dt) + 1;
t = [0 : nt-1] * dt
if nargin < 7
trace = 0;
end
if trace
trace = max((nt - 1) / 32, 1);
end
% Initialize solution, and set initial data ...
psi = zeros(nt, nx);
psire = zeros (nt, nx);
psiim = zeros (nt, nx);
psimod = zeros (nt, nx);
prob = zeros (nt, nx);
if idtype == 0 %exact family
psi(1, :) = sin(idpar * pi * x);
%elseif idtype == 1 %boosted gaussian
%idpar(1) = x0;
%idpar(2) = delta
% psi(1, :) = exp(1i*p*x)*exp(-((x - x0) ./ delta) .^ 2);
else
fprintf('diff_1d_cn: Invalid idtype=%d\n', idtype);
return
end
% Initialize Potential data
v = zeros(nx);
if vtype == 0 % No Potential
v = zeros(nx);
vpar = [0];
% elseif vtype == 1 % square barrier or well
% vpar(1) = xmin;
% vpar(2) = xmas;
% vpar(3) = vc;
% v = zeros(nx) + vc;
else
fprintf('diff_1d_cn: Invalid vtype=%d\n', vtype);
return
end
% Initialize storage for sparse matrix and RHS ...
dl = zeros(nx,1);
d = zeros(nx,1);
du = zeros(nx,1);
f = zeros(nx,1);
% Set up tridiagonal system ...
dl = ((-1i * dt) / (2 * dx^2)) * ones(nx, 1);
d = (1 + ((1i * dt) / 2) * ( v + 2.0 / dx^2)) * ones(nx,1);
du = dl;
dl2 = ((1i * dt) / (2 * dx^2)) * ones(nx, 1);
d2 = (1 - ((1i * dt) / 2) * ( v + 2.0 / dx^2)) * ones(nx,1);
du2 = dl2;
% Fix up boundary cases ...
d(1) = 1.0;
du(2) = 0.0;
dl(nx-1) = 0.0;
d(nx) = 1.0;
% Define sparse matrix ...
A = spdiags([dl d du], -1:1, nx, nx); % LHS
%B = spdiags([dl2 d2 du2], -1:1, nx, nx); % RHS
%r = zeros(nx,1);
% Compute solution using CN scheme ...
for n = 1 : nt-1
% Define RHS of linear system ...
f(2:nx-1) = psi(n, 2:nx-1) + (1i * dt * 0.5) * (( ...
psi(n, 1:nx-2) - 2 * psi(n, 2:nx-1) + psi(n, 3:nx)) / dx^2 - psi(n, 2:nx-1)+ v);
f(1) = 0.0;
f(nx) = 0.0;
% Solve system, thus updating approximation to next time
% step ...
psi(n+1, :) = A \ f;
if trace && ~mod(n,trace)
fprintf('diff_1d_cn: Step %d of %d\n', n, nt);
end
end
end
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
When I tried to implement a script to run this function:
% tdiff_1d_cn.m: Driver for diff_1d_cn ... Solution of 1d diffusion
% equation using O(dt^2, dx^2) implicit Crank-Nicolson scheme.
more off;
format long;
% idtype = 0 -> sine initial data (exact solution known)
% idtype = 1 -> gaussian initial data
% vtype = 0 : No potential
% vtype = 1 : Square welll
idtype = 0;
idpar = [1];
vtype = 0;
vpar = 0;
m = 1;
tmax = 0.2;
x0 = 0.5;
delta = 0.05;
omega = 2 * pi;
minlevel = 6;
maxlevel = 9;
olevel = 6;
ofreq = 1;
lambda = 0.1;
% Enable for MATLAB surface plots.
plotit = 1;
if plotit
close all
end
% Perform computation at various levels of discretization, store
% results in cell arrays ...
for l = minlevel : maxlevel
tstart = tic;
% Compute the solution ...
[x{1}, t{1}, psi{1}, psire{1}, psiim{1}, psimod{1}, prob{1}, v{1}] = sch_1d_cn(tmax,...
l, lambda, idtype, idpar, vtype, vpar);
telapsed = toc(tstart)
[nt{l}, nx{l}] = size(psi{l});
stride{l} = ofreq * 2^(l - olevel);
% If possible, compute exact solution ...
if idtype == 0
psixct{l} = zeros(nt{l}, nx{l});
for n = 1 : nt{l}
psixct{l}(n,:) = sin(idpar * pi * x{l});
end
end
if plotit
for it = 1 : nt{l}
figure(l);
plot(x{l},psi{l}(it,:));
xlim([0, 1]);
ylim([-1, 1]);
drawnow;
end
figure(l+10);
surf(x{l}, t{l}(1:stride{l}:end), u{l}(1:stride{l}:end, :));
drawnow;
end
% If possible, compute and output error norm ...
if idtype == 0
psierr = psi{l} - psixct{l};
fprintf('Level: %d ||error||: %25.16e\n', ...
l, sqrt(sum(sum((psierr .* psierr))) / prod(size(psierr))));
end
end
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
I got this Error:
Unable to perform assignment because
brace indexing is not supported for
variables of this type.
Error in tsch_1d_cn (line 41)
[x{1}, t{1}, psi{1}, psire{1},
psiim{1}, psimod{1}, prob{1}, v{1}]
= sch_1d_cn(tmax,...

Answers (1)

Kaashyap Pappu
Kaashyap Pappu on 24 Jan 2020
Edited: Kaashyap Pappu on 24 Jan 2020
The variables being initialized in this section of code:
psi = zeros(nt, nx);
psire = zeros (nt, nx);
psiim = zeros (nt, nx);
psimod = zeros (nt, nx);
prob = zeros (nt, nx);
are all matrices. You can access these variables using the regular brackets.
You can avoid the error by using:
[x(1), t(1), psi(1), psire(1),psiim(1), psimod(1), prob(1), v(1)] = sch_1d_cn(tmax,...
Brace indexing is used for Cell Arrays and not for regular matrices, which is why the error is being popped up.
More information about Cell Arrays can be found here: https://in.mathworks.com/help/matlab/cell-arrays.html
Hope this helps!

Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!