78
T. Torsvik
c = 1.0;
% wave speed
h = L/(M-1);
% grid-spacing
k = 0.05;
% timestep increment
for i=1:M
% Initial condition
if (i * h < 1.0)
f(i) = 0.5 * (1.0 - cos(2 * pi * i * h));
else
f(i) = 0.0;
end
end
x = 0.0:h:L;
% x-vector
plot(x,f);
% plot initial state
pause(0.1);
for j=1:N
% time step
f_new(1) = 0.0;
% boundary condition:
for i=2:M
%
f(0,t) = 0.0
dfdx = (f(i) - f(i-1))/h;
% spatial derivative
f_new(i) = f(i) - k * c * dfdx;
% time derivative
end
f = f_new;
% update old f from f_new
plot(x,f);
% plot each time step
pause(0.1);
end
There are basically two parts of the script; the first part defines the parameter
values and initial conditions (the boundary conditions are defined implicitly as zero
values), the second part solves the finite difference scheme by iterative steps in time
and space. Since the stencil of the method we use only require values at the (n) and
(n + 1) time levels we only need two arrays to define the function values for all grid
points; f stores the known function values at the (n) time step level, and f_new
stores the calculated values for the (n + 1) time step level. After each completed
sweep updating the f_new array, the old f array is replaced by the updated values.
Figure 3.9 shows the results of the numerical simulations for different values of
the time step t. These results are remarkably different from each other. In this
case we know what the result should be from the analytical solution of the PDE.
The initial profile defined by f 0 should simply be transported along characteristics
without any change in shape, which is what we see in Fig. 3.9b. Figure 3.9a shows a
transport process, but the initial profile is being smoothed out with time, which is an
effect of numerical diffusion. Figure 3.9c also shows a transport process initially, but
at the end of the simulation the picture is dominated by oscillations on the length
scale of the spatial grid resolution x which grow in amplitude with increasing