I recently asked a question about this topic, but too much was unclear to me at the time which yielded poor response. (You can find the previous question here). Anyway, I am now trying to solve the linear scalar advection equation analytically with solid wall boundary conditions (i.e. reflection at the boundaries). Something does not quite work in my code: I have a wave moving leftwards toward the left boundary, but suddenly the wave just diminshes by some factor. Then the wave proceeds to the boundary and is reflected as should be. See the results below.
Now I know what causes this problem actually, but I do not now how to resolve it. In my solution I let everything on the blue line travel to the left (propagate to the left) and then reflect on the left boundary. Thus there are waves of -0.5 amplitude going to the right "colliding" with waves going to the left. I solve this by just superpositioning these (two) waves and dividing by two. This works fine for when the colliding waves have the same amplitude, but when they are not of the same amplitude the total wave diminishes like in the picture. But how do I go about conserving amplitude or area with all these reflecting waves? What am I doing wrong?
You may test my code in MATLAB:
L = 10; H = 1; g = 9.61; w = 0.4; eps = 0.4; lambda1 = -sqrt(g*h0(L,H,w,eps)); lambda2 = sqrt(g*h0(L,H,w,eps)); N = 100; dx = L/N; xgrid = 0:dx:L; dt = 0.01; time = 0:dt:2; T = L/(-2*lambda1); for t = time(1:100) sol = zeros(1,N+1); i = 0; if 0 <= t && t < T for x = xgrid i = i + 1; if -lambda1*t <= x && x <= L + lambda1*t sol(i) = z10(x-lambda1*t,L,H,w,eps); elseif 0 <= x && x < -lambda1*t xarg = x-lambda1*t; sol(i) = (1/2)*(zr(x,t,L,H,w,eps,lambda1) + z10(xarg,L,H,w,eps)); elseif L+lambda1*t < x && x <= L sol(i) = z10(L,L,H,w,eps); end end elseif T <= t && t <= 2*T for x = xgrid i = i + 1; if L+lambda1*t <= x && x <= -lambda1*t sol(i) = zr(x,t,L,H,w,eps,lambda1); elseif 0 <= x && x < L+lambda1*t sol(i) = (1/2)*(zr(x,t,L,H,w,eps,lambda1) + z10(x-lambda1*t,L,H,w,eps)); elseif -lambda1*t <= x && x <= L sol(i) = z10(L,L,H,w,eps); end end end plot(xgrid,sol) title(['t=',num2str(t)]) axis([0 L -0.7 -0.5]) pause(0.1) end plot(xgrid,sol) xlabel('x') ylabel('q') axis([0 L -0.7 -0.5]) title(['Analytical solution for linearized system at t=',num2str(t)]) function val = h0(L,H,w,eps) val = H + (eps*w*sqrt(pi))*erf(L/(2*w))/L; end function val = zr(x,t,L,H,w,eps,lambda1) % Right boundary at x=0 tr = t + x/lambda1; xarg = -lambda1*tr; val = z10(xarg,L,H,w,eps); end function val = z10(x,L,H,w,eps) % Initial data val = -(H + eps*exp(-(x-L/2)^2/(w^2)))/2; end
The code "works" for t < 3, but I have not implemented boundary conditions at the other end so it will not show correctly for times larger than that. I hope I make the question easy to answer!