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 \]
Using initial data \[ u_0(x) = \sin(2\pi x) \] the exact solution is \[ u(x,t) = e^{-4\pi^2t}\sin(2\pi x) \]
The Crank-Nicolson FDA is
\[ \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 \]
This can be written in the form \[ c_+ u^{n+1}_{j+1} + c_0 u^{n+1}_j + c_- u^{n+1}_{j-1} = f_j \] where \[ c_+ = -\frac{1}{2\Delta x^2} \] \[ c_0 = \frac{1}{\Delta t} + \frac{1}{\Delta x^2} \] \[ c_- = -\frac{1}{2\Delta x^2} = c_+ \]
My implementation of the solution as the function diff_1d_cn is as
follows:
function [x t u] = diff_1d_cn(tmax, level, dtbydx, ...
omega, x0, delta, idtype, trace)
% diff_1d_imp: Solves 1D diffusion equation using O(dt^2,dx^2) Crank-Nicolson
% implicit scheme. Diffusion constant, sigma = 1.
%
% Inputs
%
% tmax: Maximum integration time.
% level: Spatial discretization level.
% dtbydx: Ratio of spatial and temporal mesh spacings.
% omega: Initial data parameter (sinusoidal data).
% x0: Initial data parameter (Gaussian data).
% delta: Initial data parameter (Gaussian data).
% idtype: Selects initial data type.
% 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 = dtbydx * dx;
nt = round(tmax / dt) + 1;
t = [0 : nt-1] * dt;
if nargin < 8
trace = 0;
end
if trace
trace = max((nt - 1) / 32, 1);
end
% Initialize solution, and set initial data ...
u = zeros(nt, nx);
if idtype == 0
u(1, :) = sin(omega * x);
elseif idtype == 1
u(1, :) = exp(-((x - x0) ./ delta) .^ 2);
else
fprintf('diff_1d_cn: Invalid idtype=%d\n', idtype);
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 = -0.5 / dx^2 * ones(nx, 1);
d = (1.0 / dt + 1.0 / dx^2) * ones(nx,1);
du = dl;
% 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);
% Compute solution using CN scheme ...
for n = 1 : nt-1
% Define RHS of linear system ...
f(2:nx-1) = u(n, 2:nx-1) / dt + 0.5 * ( ...
u(n, 1:nx-2) - 2 * u(n, 2:nx-1) + u(n, 3:nx)) / dx^2;
f(1) = 0.0;
f(nx) = 0.0;
% Solve system, thus updating approximation to next time
% step ...
u(n+1, :) = A \ f;
if trace && ~mod(n,trace)
fprintf('diff_1d_cn: Step %d of %d\n', n, nt);
end
end
end
Here is a script tdiff_1d_cn that, among other things, computes the
solution error norms as a function of discretization level.
% 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
idtype = 0
tmax = 0.2
x0 = 0.5
delta = 0.05
omega = 2 * pi
minlevel = 6
maxlevel = 9
olevel = 6
ofreq = 1
dtbydx = 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{l}, t{l}, u{l}] = diff_1d_cn(tmax, l, dtbydx, omega, ...
x0, delta, idtype, 0);
telapsed = toc(tstart)
[nt{l}, nx{l}] = size(u{l});
stride{l} = ofreq * 2^(l - olevel);
% If possible, compute exact solution ...
if idtype == 0
uxct{l} = zeros(nt{l}, nx{l});
for n = 1 : nt{l}
uxct{l}(n,:) = exp(-omega^2 * t{l}(n)) * sin(omega * x{l});
end
end
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
% If possible, compute and output error norm ...
if idtype == 0
uerr = u{l} - uxct{l};
fprintf('Level: %d ||error||: %25.16e\n', ...
l, sqrt(sum(sum((uerr .* uerr))) / prod(size(uerr))));
end
end
Partial output from the script is:
Level: 6 ||error||: 6.0516904501925616e-05
telapsed =
0.012512000000000
Level: 7 ||error||: 1.5205493653807096e-05
telapsed =
0.019560000000000
Level: 8 ||error||: 3.8117334990173340e-06
telapsed =
0.056418000000000
Level: 9 ||error||: 9.5428113131619461e-07
from which it can be seen that the error norms are going down by very nearly a factor of 4 for each halving of the mesh spacings, indicating \( O(\Delta t^2, \Delta x^2) = O(h^2) \) convergence.
We now want to consider
\[ 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 we discussed in class is (in the form appropriate for numerical solution): \[ 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} \] Let us consider 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 initial 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) \] \[ u^2_j = u^1_j + \frac{1}{2} \Delta t^2 (u''_0)_j \] When implementing the last expression, you should use the usual \( O(h^2) \) approximation of the second derivative to evaluate \( u''_0 \) rather than computing it analytically (simply since it's arguably more straightforward and sufficiently accurate).
EXERCISE
Write a MATLAB function, wave_1d, with header and function arguments
as follows
function [x t u] = wave_1d(tmax, level, lambda, x0, delta)
% 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).
%
% Outputs
%
% x: Discrete spatial coordinates.
% t: Discrete temporal coordinates.
% u: Computed solution.
wave_1d is to solve the FDA described above. You may find it useful
to use the code for diff_1d_cn as a template.
DOWNLOAD SOURCE
In addition to wave_1d, write a driver script twave_1d which calls
wave_1d with appropriate parameters and makes plots of the solution
(loop over the time coordinate) if plotting is enabled.
To start, I suggest that you experiment with the following parameter settings:
tmax = 3.0level = 6 through 11lambda = 0.5x0 = 0.4delta = 0.05Now try a run with
tmax = 3.0level = 8lambda = 1.01x0 = 0.4delta = 0.05What do you observe?
If you have time, convergence test your implementation by augmenting your script (or writing a new one) so that it
wave_1d and saving the results.What can you say about the convergence of your solutions?
You may find it useful to refer to the script tdiff_1d_cn for ideas
about how to implement the driver script.
DOWNLOAD SOURCE