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.