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,...
