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.