# 2D wave equation with Mur boundary condition - setting up the RHS and solving (time-steps)

I am trying to solve a 2D wave equation implicitly using FD with central approximations with the following boundary conditions

\begin{align} &u=2\sin\left(\frac{2\pi}{5}t\right)\quad \text{at }x=0\\ &\frac{\partial u}{\partial y}=0 \quad \text{at }y=0\\ &\frac{\partial u}{\partial y}=0 \quad \text{at }y=2\\ &\frac{\partial u}{\partial t}=-\frac{\partial u}{\partial x} \quad \text{at }x=5 \end{align}

graphically represented as The boundary condition at x=0 generates the wave. The boundary condition at x=5 refers to Mur boundary condition, i.e. the wave doesn't bounce back, but simply continues to move outside the domain. The general equation is given as The Mur boundary condition can be expressed mathematically as The boundary conditions were evaluated analytically and these are the 7 equations following the entire domain: Problem- I have managed to write some code, including setting up the [A] matrix (the LHS of the equations), but I have problems with formulating the [b] matrix (RHS) and solving it through the time-steps. I would appreciate if someone show me how to write it for the first equation - then I will be able to do it dor the remaining 6.

% length, time, height
L = 5; % [m]
h = 2; % [m]

d = 0.25; % space spacing
dt = 0.05; % time increment
M = 100;
T_max = M*dt; % [s]

nx = floor(L/d); % number of <x> samples
ny = floor(h/d); % number of <y> samples
k = floor(T_max/dt) + 1; % number of time samples

% Constants
alpha = dt/d;
r = dt^2/d^2;

% Number of grid points:
N=nx*ny;

% Initialise matrices
U = zeros(N,1);
b = zeros(N,1);
x = zeros(N,1);
A = zeros(N,N);

% Set numbering for the [A] matrix
num = 1;
for i=1:nx
for j=1:ny
number(i,j) = num;
num = num + 1;
end
end

% [A] matrix - boundary conditions
% Case 1
for i=2:nx-1
for j=2:ny-1
ii=number(i,j);
A(ii, number(i, j)) = 1+4*r;
A(ii, number(i+1, j)) = -r;
A(ii, number(i-1,j)) = -r;
A(ii, number(i, j+1)) = -r;
A(ii, number(i, j-1)) = -r;
end
end

% Case 2
for i=1
for j=2:ny-1
ii=number(i,j);
A(ii, number(1, j)) = 1+4*r;
A(ii, number(2, j)) = -r;
A(ii, number(1, j+1)) = -r;
A(ii, number(1, j-1)) = -r;
end
end

% Case 3
for i=nx
for j=1:ny
ii=number(i,j);
A(ii, number(nx, j)) = 1+alpha;
A(ii, number(nx-1,j)) = 1-alpha;
end
end

% Case 4
for i=2:nx-1
for j=1
ii=number(i,j);
A(ii, number(i, j)) = 1+4*r;
A(ii, number(i-1,j)) = -r;
A(ii, number(i+1,j)) = -r;
A(ii, number(i, j+1)) = -2*r;
end
end

% Case 5
for i=1
for j=1
ii=number(i,j);
A(ii, number(i,j)) = 1+4*r;
A(ii, number(i+1,j)) = -r;
A(ii, number(i, j+1)) = -2*r;
end
end

% Case 6
for i=2:nx-1
for j=ny
ii=number(i,j);
A(ii, number(i,j)) = 1+4*r;
A(ii, number(i-1,j)) = -r;
A(ii, number(i+1, j)) = -r;
A(ii, number(i, j-1)) = -2*r;
end
end

% Case 7
for i=1
for j=ny
ii=number(i,j);
A(ii, number(i,j)) = 1+4*r;
A(ii, number(i+1, j)) = -r;
A(ii, number(i, j-1)) = -2*r;
end
end

% Inverse of [A] matrix
A_inv = inv(A);

% [b] matrix (RHS)
time = 0;
for t=1:M
% Computing [b] matrix
for i=1:nx
for j=1:ny
ii = number(i,j);
% U(1,1) = 2*sin((2*pi/5)*(t*dt));

% Case 1
for i=2:nx-1
for j=2:ny-1

end
end
% or

for i=1:nx
for j=1:ny

end
end

end
end

% The solution vector
x = A_inv*b;

% 2D matrix from the solution vector

• Welcome to SciComp Exchange! Can you provide the (non discrete) formulation of the Mur BCs? Also, can you use MathJax for the equations? (I already formatted the first set of them). – nicoguaro Apr 22 '15 at 19:49