10.5 Exercises 327
doublealpha = tstep*tstep/(h*h)
// We define the solution u at an explicit time step l
// using Armadillotodefine matrices
mat u( n+1, n+1), u_last( n+1, n+1), u_next( n+1, n+1);
// We declare also vectors that hold theposition inthe x and y directions
vec x(n+1), y(n+1);
for( int i = 0; i < n+1 ; i++ ){
x(i) = i*h;
y(i) = x(i);
}
u_last = 0.0;
// initializing thefunction(using the initial conditions)
for( int i = 1; i < n; i++ ){// setting initial step
for( int j = 1; j < n; j++ ){
u_last(i,j) = sin(PI*x(i))*sin(2*PI*y(j));
}
}
u = 0.0; u_next =0.0; // includes also the boundary, setforall times
for( int i = 1; i < n; i++ ){// setting first step using the initial derivative
for( int j = 1; j < n; j++ ){
u(i,j) = u_last(i,j) - alpha*0.5*
(4*u_last(i,j) - u_last(i+1,j) - u_last(i-1,j) - u_last(i,j+1) - u_last(i,j-1));
}
// iteratingintime
doublet = tinitial;
while( t < tfinal ){
t = t + tstep;
for( int i = 1; i < n; i++ ){// computing next step
for( int j = 1; j < n ; j++ ){
u_next(i,j) = 2*u(i,j) - u_last(i,j) -
alpha*(4*u(i,j) - u(i+1,j) - u(i-1,j) - u(i,j+1) - u(i,j-1));
}
}
// Updatethenthe newfunctionvalue
for( int i = 1; i < n; i++ ){
for( int j = 1; j < n; j++ ){
u_last(i,j) = u(i,j);
u(i,j) = u_next(i,j);
}
}
// One may considerto print to filethe results after a given # of time steps
.... printstatements
}
return0;
}We modify now the wave equation in order to consider a 2 + 1 dimensional wave equation
with a position dependent velocity, given by
∂^2 u
∂t^2
=∇·(λ(x,y)∇u).Ifλis constant, we obtain the standard wave equation discussedin the two previous points.
The solutionu(x,y,t)could represent a model for water waves. It represents then the surface
elevation from still water. We will modelλas
λ=gH(x,y),withgbeing the acceleration of gravity andH(x,y)is the still water depth.