We are solving
\[ u_{tt} = u_{xx} \] on \[ 0 \le x \le 1 \quad\quad t \ge 0 \] subject to \[ u(x,0) = u_0(x) \] \[ u_t(x,0) = v_0(x) \] \[ u(0,t) = u(1,t) = 0 \] The FDA is \[ u^{n+1}_j = 2u^n_j - u^{n-1}_j + \lambda^2 \left( u^n_{j+1} - 2 u^n_j + u^n_{j-1}\right) \quad\quad j = 2, 3, \ldots n_x - 1 \quad\quad n = 2, 3, \ldots \] \[ u^{n+1}_1 = u^{n+1}_{n_x} = 0 \] where \[ \lambda = \frac{\Delta t}{\Delta x} \] We are considering the specific initial data \[ u_0(x; x_0, \delta) = e^{-((x-x_0)/\delta)^2} \] \[ v_0(x) = 0 \] which is so-called time symmetric data due to the vanishing of \(v_0\).
The initialization of the FDA is then \[ u^1_j = u_0(x_j; x_0,\delta) \] and using Taylor series expansion \[ u^2_j = u^1_j + \Delta t \left( u_t \right)^1_j + \frac{1}{2} \Delta t^2 (u_{tt})^1_j + O(\Delta t^3) \] Now, \( (u_t)^1_j=0 \) since \( v_0 = 0 \), and \( u_{tt} = u_{xx} \) from the PDE. Thus, to \( O(\Delta t^2) \) accuracy we have \[ u^2_j = u^1_j + \frac{1}{2} \Delta t^2 (u''_0)_j \]
wave_1dfunction [x t u] = wave_1d(tmax, level, lambda, x0, delta, trace)
% wave_1d: Solves 1d wave equation using O(dt^2,dx^2) explicit scheme.
%
% Inputs
%
% tmax: Maximum integration time.
% level: Spatial discretization level.
% lambda: Ratio of spatial and temporal mesh spacings.
% x0: Initial data parameter (Gaussian data).
% delta: Initial data parameter (Gaussian data).
% trace: Controls tracing output.
%
% Outputs
%
% x: Discrete spatial coordinates.
% t: Discrete temporal coordinates.
% u: Computed solution.
% Define mesh and derived parameters ...
nx = 2^level + 1;
x = linspace(0.0, 1.0, nx);
dx = x(2) - x(1);
dt = lambda * dx;
nt = round(tmax / dt) + 1;
t = [0 : nt-1] * dt;
if nargin < 6
trace = 0;
end
if trace
trace = max((nt - 1) / 32, 1);
end
% Initialize solution matrix, and set initial data ...
% Initial data is time symmetric (u_t(x,0) = 0) Gaussian profile ...
u = zeros(nt, nx);
u(1, 2:nx-1) = exp(-((x(2:nx-1) - x0) ./ delta) .^ 2);
% Second time step comes from Taylor expansion, PDE (u_tt = u_xx),
% and numerical approximation to second derivative ...
u(2, 2:nx-1) = u(1, 2:nx-1) + ...
0.5 * lambda^2 * (u(1, 1:nx-2) - 2.0 * u(1, 2:nx-1) + u(1, 3:nx));
% Compute solution using FDA ...
for n = 2 : nt-1
u(n+1, 2:nx-1) = 2.0 * u(n, 2:nx-1) - u(n-1, 2:nx-1) + ...
lambda^2 * (u(n, 1:nx-2) - 2.0 * u(n, 2:nx-1) + u(n, 3:nx));
if trace && ~mod(n,trace)
fprintf('diff_1d_cn: Step %d of %d\n', n, nt);
end
end
end
twave_1d% twave_1d.m: Driver for wave_1d ... Solution of 1d wave
% equation using O(dt^2, dx^2) explicit scheme (standard scheme) ...
more off;
format long;
tmax = 3.0
x0 = 0.4
delta = 0.05
minlevel = 6
maxlevel = 10
olevel = 6
ofreq = 1
lambda = 0.5
% Enable for plotting ...
plotit = 0
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{l}, t{l}, u{l}] = wave_1d(tmax, l, lambda, x0, delta, 1);
telapsed = toc(tstart)
[nt{l}, nx{l}] = size(u{l});
stride{l} = ofreq * 2^(l - olevel);
if plotit
for it = 1 : nt{l}
figure(l);
plot(x{l},u{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
end
% Perform convergence analysis ...
if maxlevel - minlevel > 2
for l = minlevel : maxlevel - 1
du = u{l+1}(1:2:end, 1:2:end) - u{l};
fprintf('||u(%d)-u(%d)||: %g\n', l+1, l, norm(du)/sqrt(nx{l} * nt{l}));
end
end
twave_1d||u(7)-u(6)||: 0.0205912
||u(8)-u(7)||: 0.00544933
||u(9)-u(8)||: 0.00136913
||u(10)-u(9)||: 0.000342646
From these results we see that the solution is converging as an \( O(h^2) \) quantity, as expected.
When \( \lambda > 1 \) we expect---on the basis of the von Neumann analysis performed in class---that our explicit FDA will become unstable, and this is indeed what is seen. Specifically, the solution develops grid point to grid point oscillations which grow extremely rapidly.
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.
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?