Finite difference method solves the wave equation for one set of bounday conditions, but does not when I change them

33 Views Asked by At

I wrote a finite difference algorithm in Matlab to solve the wave equation which is derived here.

When I ran my code, the plotted graphs of the numerical and analytical solution deviated, which is the problem I am trying to solve. The finite difference algorithm is given in the snippet.

t_min = 0;
t_max = 10;
rps   = 400;                 % Resolution per second
nt    = t_max * rps;
t     = linspace(t_min, t_max, nt);
dt    = t(2) - t(1);

x_min = 0;
x_max = 8;
rpm   = 100;                  % Resolution per meter
nx    = x_max * rpm;
x     = linspace(x_min, x_max, nx); 
dx    = x(2) - x(1);

c     = 3;                    % Wave speed
A     = pi;                   % Amplitude        
L     = x_max;                % Rod lenght   
w_o   = 1;                    % Circular frequency

Cn = c * (dt / dx);           % Courrant number

if Cn > 1                     % Stability criteria

    error('The stability condition is not satisfied.');
end

U = zeros(nx, nt);           

U(:, 1)   = zeros(nx, 1);    
U(:, 2)   = zeros(nx, 1);     

U(1, :)   = A * sin(w_o * t);
U(end, :) = zeros(1, nt);    

for j = 2 : (nt - 1)

    for i = 2 : (nx - 1)
 
        U(i, j + 1) = Cn^2 * ( U(i + 1, j) - 2 * U(i, j) + U(i - 1, j) ) + 2 * U(i, j) - U(i, j - 1);
    end     
end

figure('Name', 'Numeric solution');
surface(t, x, U, 'edgecolor', 'none');
grid on
colormap(jet(256));
colorbar;
xlabel('t ({\its})');
ylabel('x ({\itm})');
title('U(t, x) ({\itm})');

To find the bug, I tried to change the boundary conditions and see if my graph would look better. It turned out that it did, which means that the code in my double for loop is ok. The boundary conditions are the problem. However, I do not know why the code works with the new boundary conditions and does not with the old ones. I am hoping that somebody will point this out to me. The code I ran is given in the snippet.

t_min = 0;
t_max = 1;
rps   = 400;                 % Resolution per second
nt    = t_max * rps;
t     = linspace(t_min, t_max, nt);
dt    = t(2) - t(1);

x_min = 0;
x_max = 1;
rpm   = 100;                  % Resolution per meter
nx    = x_max * rpm;
x     = linspace(x_min, x_max, nx); 
dx    = x(2) - x(1);

c     = 3;                    % Wave speed
A     = pi;                   % Amplitude        
L     = x_max;                % Rod lenght   
w_o   = 1;                    % Circular frequency

Cn = c * (dt / dx);           % Courrant number

if Cn > 1                     % Stability criteria

    error('The stability condition is not satisfied.');
end

U = zeros(nx, nt);           

U(:, 1)   = sin(pi*x);     
U(:, 2)   = sin(pi*x) * (1 + dt);   

U(1, :)   = zeros(1, nt);
U(end, :) = zeros(1, nt);  

for j = 2 : (nt - 1)

    for i = 2 : (nx - 1)
 
        U(i, j + 1) = Cn^2 * ( U(i + 1, j) - 2 * U(i, j) + U(i - 1, j) ) + 2 * U(i, j) - U(i, j - 1);
    end
 end

figure('Name', 'Numeric solution');
surface(t, x, U, 'edgecolor', 'none');
grid on
colormap(jet(256));
colorbar;
xlabel('t ({\its})');
ylabel('x ({\itm})');
title('U(t, x) ({\itm})');

To summarize, I am trying to understand why solving the wave equation with the boundary conditions:

$$ U(0,x)=0 $$ $$ U_t(0,x)=0 $$ $$ U(t,0)=A \sin(w_0t) $$ $$ U(t,x_{max})=0 $$

gives incorrect results when using finite differences, but solving the wave equation with the boundary conditions:

$$ U(0,x)=\sin (\pi x) $$ $$ U_t(0,x)=\sin (\pi x) $$ $$ U(t,0)=0 $$ $$ U(t,x_{max})=0 $$

gives correct results. Can anybody explain this to me?

I asked this question on stack overflow, but got redirected here.