function u2 = fd1d_wave_start ( x_num, x_vec, t, t_delta, alpha, u_x1, u_x2, ...
ut_t1, u1 )
%*****************************************************************************80
%
%% fd1d_wave_start() takes the first step for the wave equation.
%
% Discussion:
%
% This program solves the 1D wave equation of the form:
%
% Utt = c^2 Uxx
%
% over the spatial interval [X1,X2] and time interval [T1,T2],
% with initial conditions:
%
% U(T1,X) = U_T1(X),
% Ut(T1,X) = UT_T1(X),
%
% and boundary conditions of Dirichlet type:
%
% U(T,X1) = U_X1(T),
% U(T,X2) = U_X2(T).
%
% The value C represents the propagation speed of waves.
%
% The program uses the finite difference method, and marches
% forward in time, solving for all the values of U at the next
% time step by using the values known at the previous two time steps.
%
% Central differences may be used to approximate both the time
% and space derivatives in the original differential equation.
%
% Thus, assuming we have available the approximated values of U
% at the current and previous times, we may write a discretized
% version of the wave equation as follows:
%
% Uxx(T,X) = ( U(T, X+dX) - 2 U(T,X) + U(T, X-dX) ) / dX^2
% Utt(T,X) = ( U(T+dt,X ) - 2 U(T,X) + U(T-dt,X ) ) / dT^2
%
% If we multiply the first term by C^2 and solve for the single
% unknown value U(T+dt,X), we have:
%
% U(T+dT,X) = ( C^2 * dT^2 / dX^2 ) * U(T, X+dX)
% + 2 * ( 1 - C^2 * dT^2 / dX^2 ) * U(T, X )
% + ( C^2 * dT^2 / dX^2 ) * U(T, X-dX)
% - U(T-dT,X )
%
% (Equation to advance from time T to time T+dT, except for FIRST step!)
%
% However, on the very first step, we only have the values of U
% for the initial time, but not for the previous time step.
% In that case, we use the initial condition information for dUdT
% which can be approximated by a central difference that involves
% U(T+dT,X) and U(T-dT,X):
%
% dU/dT(T,X) = ( U(T+dT,X) - U(T-dT,X) ) / ( 2 * dT )
%
% and so we can estimate U(T-dT,X) as
%
% U(T-dT,X) = U(T+dT,X) - 2 * dT * dU/dT(T,X)
%
% If we replace the "missing" value of U(T-dT,X) by the known values
% on the right hand side, we now have U(T+dT,X) on both sides of the
% equation, so we have to rearrange to get the formula we use
% for just the first time step:
%
% U(T+dT,X) = 1/2 * ( C^2 * dT^2 / dX^2 ) * U(T, X+dX)
% + ( 1 - C^2 * dT^2 / dX^2 ) * U(T, X )
% + 1/2 * ( C^2 * dT^2 / dX^2 ) * U(T, X-dX)
% + dT * dU/dT(T, X )
%
% (Equation to advance from time T to time T+dT for FIRST step.)
%
% It should be clear now that the quantity ALPHA = C * DT / DX will affect
% the stability of the calculation. If it is greater than 1, then
% the middle coefficient 1-C^2 DT^2 / DX^2 is negative, and the
% sum of the magnitudes of the three coefficients becomes unbounded.
%
% Licensing:
%
% This code is distributed under the MIT license.
%
% Modified:
%
% 24 January 2012
%
% Author:
%
% John Burkardt
%
% Reference:
%
% George Lindfield, John Penny,
% Numerical Methods Using MATLAB,
% Second Edition,
% Prentice Hall, 1999,
% ISBN: 0-13-012641-1,
% LC: QA297.P45.
%
% Input:
%
% integer X_NUM, the number of nodes in the X direction.
%
% real X_VEC(X_NUM), the spatial coordinates of the nodes.
%
% real T, the time after the first step has been taken.
% In other words, T = T1 + T_DELTA.
%
% real T_DELTA, the time step.
%
% real ALPHA, the stability coefficient, computed by FD1D_WAVE_ALPHA.
%
% real U_X1(T), U_X2(T), functions for the left and right boundary
% conditions.
%
% real UT_T1(X), the function that evaluates dUdT at the initial time.
%
% real U1(X_NUM), the initial condition.
%
% Output:
%
% real U2(X_NUM), the solution at the first time step.
%
ut = ut_t1 ( x_num, x_vec );
u2 = zeros ( 1, x_num );
u2(1) = u_x1 ( t );
u2(2:x_num-1) = alpha^2 * u1(3:x_num) / 2 ...
+ ( 1 - alpha^2 ) * u1(2:x_num-1) ...
+ alpha^2 * u1(1:x_num-2) / 2 ...
+ t_delta * ut(2:x_num-1);
u2(x_num) = u_x2 ( t );
return
end