fcn_deutn.mfcn_deutn incorporates the ODE for normalization of the
wave function.
function dydx = fcn_deutn(x, y)
% Function fcn_deut evaluates derivatives for toy deuteron problem,
% including normalization.
global x0 E;
dydx = ones(3,1);
dydx(1,1) = y(2);
if x <= x0
dydx(2,1) = (-1 - E) * y(1);
else
dydx(2,1) = -E * y(1);
end
dydx(3,1) = y(1)^2;
end
deutn_bisect.m% deutn_bisect: Solves BVP for toy deuteron problem with wavefunction
% normalization and bisection
global x0 E;
% Domain outer boundary ...
xmax = 60.0;
% Tolerance parameters for ode45 ...
abstol = 1.0e-8;
reltol = 1.0e-8;
options = odeset('AbsTol', abstol * [1 1 1]', 'RelTol', reltol);
% Tolerance parameter for bisection search
Etol = 1.0e-10;
% Parameters ...
x0 = 2.0;
% Initial bracket ...
Emin = -0.5;
Emax = -0.01;
E = Emin;
[xout yout] = ode45(@fcn_deutn, [0.0 xmax], [0.0 1.0 0.0]', options);
fmin = yout(end,1);
% Perform bisection search on eigenvalue ...
epsilon = 2.0 * Etol;
while epsilon > Etol
E = 0.5 * (Emin + Emax)
[xout yout] = ode45(@fcn_deutn, [0.0 xmax], [0.0 1.0 0.0]', options);
clf;
plot(xout, yout(:,1));
drawnow;
xx = input('Continue? ');
fmid = yout(end,1);
if fmid == 0
break;
elseif fmid * fmin < 0
Emax = E
else
Emin = E
fmin = fmid;
end
if E == 0
epsilon = abs(Emax - Emin);
else
epsilon = abs((Emax - Emin) / E);
end
end
% Normalize wave function ...
yout(:,1) = yout(:,1) / sqrt(yout(round(end/2),3));
% Check normalization ...
Emin
Emax
norm = trapz(xout(1:round(end/2)), yout(1:round(end/2),1).^2)
% Make final plot and save as JPEG ...
figure(1);
clf;
hold on;
axis square;
box on;
xlabel('x');
ylabel('u');
plot(xout(round(1:end/2)), yout(round(1:end/2), 1));
title_str = sprintf('Toy deuteron problem: x_0 = %g', x0);
title(title_str);
legend_str = sprintf('E = %g', E);
legend(legend_str);
fname = sprintf('deutn-%g.jpg', x0);
print(fname,'-djpeg');
Recall that we are solving
\[ u_t(x,t) = u_{xx} \] (unit diffusion constant) on \[ 0 \le x \le 1 \quad\quad t \ge 0 \] subject to \[ u(x,0) = u_0(x) \] and \[ u(0,t) = u(1,t) = 0 \] The scripts that you will use below are configured so that the initial data is \[ u_0(x) = \sin(2\pi x) \] corresponding to a solution \[ u(x,t) = e^{-4\pi^2t}\sin(2\pi x) \]
This scheme has \( O(\Delta t, \Delta x^2) \) truncation error.
DOWNLOAD SOURCE
EXERCISE
Run the script tdiff_1d_ftcs.m. What do you observe? Now modify
tdiff_1d_ftcs.m so that the parameter alpha is set to 0.51. Rerun
the script. What happens?
Reset alpha to 0.49, set plotit to 0 and maxlevel to 9. Rerun
the script. What can you say about the behaviour of the error norm?
Can you explain the behaviour, particularly in view of the nature of
the truncation error?
This scheme also has \( O(\Delta t, \Delta x^2) \) truncation error.
DOWNLOAD SOURCE
The driver script tdiff_1d_imp.m is configured so that as level is
varied, the ratio \(\Delta t/\Delta x\) (dtbydx) is held fixed.
EXERCISE
Run the script. What is the behaviour of the error norm in this case, and can you explain it?
Now study the function diff_id_imp and ensure that you understand how
it works, especially how the tridiagonal system is set up and solved.
Note that the linear system that needs to be solved is the
same at each time step, so we can define the sparse matrix \( A \)
outside of the time stepping loop. On the other hand, the right hand
side of the system changes from step to step (since it depends on
\( u^n \)), so its definition needs to be inside the loop.
EXERCISE
Working from diff_1d_imp.m and tdiff_1d_imp.m, implement a solution
of the diffusion equation using the Crank-Nicolson scheme:
\[ \frac{u^{n+1}_j - u^n_j}{\Delta t} = \frac 12 \left( \frac{u^{n+1}_{j+1}-2u^{n+1}_j+u^{n+1}_{j-1}}{\Delta x^2} + \frac{u^{n}_{j+1}-2u^{n}_j+u^{n}_{j-1}}{\Delta x^2} \right), \quad\quad j = 2,3,\ldots n_x -1 \] \[ u^{n+1}_1 = u^{n+1}_{n_x} = 0 \]
To get started, rewrite the above equation in the form \[ c^+_j u^{n+1}_{j+1} + c^0_j u^{n+1}_j + c^-_j u^{n+1}_{j-1} = f_j \] where \(c^+_j\), \(c^0_j\) and \(c^-_j\) are coefficients.
Note that once you've set up your tridiagonal matrix, say A,
you can see the corresponding full matrix using the MATLAB command
full(A).
What can you say about the error norm as a function of discretization level for this scheme?