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})');