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 our model 2D problem for studying relaxation methods is \[ \nabla^2 u(x,y) \equiv u_{xx} + u_{yy} = f(x,y) \] on \[ 0 \le x \le 1, \,\, 0 \le y \le 1 \] subject to \[ u(0,y) = u(1,y) = u(x,0) = u(x,1) = 0 \]
Download the script file sor2d.m and the function file
calc_exact.m.
DOWNLOAD SOURCE
The code is configured so that the exact solution is \[ u_{\rm exact}(x,y) = \sin(\omega_x x)\sin(\omega_y y) \]
where \( \omega_x \) and \( \omega_y \) are integer multiples of \(\pi\), corresponding to a source function
\[ f(x,y) = -(\omega_x^2 + \omega_y^2) u_{\rm exact}(x,y) \]
EXERCISE
Edit the script file, removing all of the code that has to do with the
generation of an AVI movie, but keep the 3 lines
uerr = uxct - u;
surf(x, y, u);
title('Sweep 0');
from before the main loop, and the three lines
surf(x, y, u);
titlestr = sprintf('Sweep %g\n', isweep);
title(titlestr);
from within the main loop, so that the script still plots the solution as the relaxation proceeds. Also add
drawnow;
after the last title command so that the graphics get displayed as
the iteration progresses.
EXERCISE
Run the script and observe its behaviour.
Now consider a modified, nonlinear problem \[ u_{xx} + u_{yy} + u^2 = f(x,y) \] on \[ 0 \le x \le 1, \,\, 0 \le y \le 1 \] subject to \[ u(0,y) = u(1,y) = u(x,0) = u(x,1) = 0 \]
EXERCISE
Modify sor2d.m and calc_exact.m to solve this problem using SOR. I
suggest that you first copy these files to new ones named sor2d_nl.m and
calc_exact_nl.m (for example).
In performing the modifications, do the following
Make sor2d_nl a function with the following header
function [u uxct] = sor2d_nl(level, rtol)
The input arguments are the discretization level, and a tolerance parameter for the residual norm. The output arguments are the computed and exact solutions, respetively.
Rather than iterating for a fixed number, nsweep, of relaxation
sweeps, have the SOR function return when the residual norm is
less that or equal to the rtol parameter.
Modify the calc_exact_nl function so that it computes and returns
the appropriate exact solution (which is the same as the original) and
right hand side (which is not the same as the original).
Modify the SOR iteration in the SOR function appropriately.
EXERCISE
Using calculations on levels 5, 6 and 7, verify that your implementation
appears to converge properly. Use rtol = 1.0e-8.
Denote by \( u^{h} \) and \( u^{2h} \) the finite difference solutions computed with mesh spacings \( h \) and \( 2h \), respectively. Equivalently we have that \( u^h \) is a level-\(l\) solution while \( u^{2h} \) is a level-\(l-1\) solution.
The FDA we are using to solve our elliptic equation is \(O(h^2)\) accurate. From this we can reasonably assume that, at least in the limit that \(h\to0\), we have
\[ \lim_{h\to0} u^h(x,y) = u(x,y) + h^2 e_2(x,y) +O(h^4) \] where \(u(x,y)\) is the continuum solution and \(e_2(x,y)\) is the leading order error function for the discrete solution.
EXERCISE
Derive a linear combination of \( u^{h} \) and \( u^{2h} \) that is equal to \( u(x,y) \) at leading order and which eliminates the \(h^2\) error term. I.e., if \(\alpha\) and \(\beta\) are judiciously chosen constants, then
\[ \alpha u^{h} + \beta u^{2h} = u(x,y) + O(h^4) \]
EXERCISE
Use your formula to extrapolate level-6 and level-7 solutions and record the error norms for those two solutions as well as that for the extrapolated solution. What can you say about the various error norms?