1

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})');
  • 1
    I think this is really not a suitable question for SO. I would ask this question somewhere like Mathematics StackExchange with more focus on the math basics of the question than the coding aspect of it. I guess you do not even need to show the code. – If_You_Say_So Nov 15 '22 at 14:47
  • I think you are right now that I think about it. Thank you for the notice. – Nikola Ristic Nov 15 '22 at 17:33

0 Answers0