Why does changing the boundary conditions cause the finite difference algorithm to diverge?

74 views Asked by At

I wrote a finite difference algorithm 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 menter
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 tryed 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 menter
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})');
0

There are 0 answers