4
$\begingroup$

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 schematics

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

general equation

The Mur boundary condition can be expressed mathematically as Mur BC

The boundary conditions were evaluated analytically and these are the 7 equations following the entire domain:

equations

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
$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – nicoguaro Apr 22 '15 at 19:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.