# Numerical instability in finite difference string simulation

I'm trying to simulate a string with closed edges (fixed at 0) using simple RK4 finite difference scheme.

Generally it seems that my simulation works as intended, except from (I think) the boundary conditions: Instead of reflecting the incoming wave perfectly, some the the wave is converted to smaller sine-like waves trailing the main wavefront. See attached video (Google Drive link).

Is there a specific name for the phenomenon? and what should I do to fix it?

Thanks a lot!

P.S. here is the simulation code in Matlab. Probably not the best code ever, but it works.

%% Aesthetics
clear all; clc; figure(1); clf;

%% Constants
C = Constants();
h = figure(1);

%% Generate space and initial state
X = repmat(linspace(0, 2, C.K), C.N, 1);

P0 = exp(-100.*(X-1).^2);
P0(2:end, :) = 0;
V0 = zeros(size(X));
V0(1,:) = 200.*(X(1,:)-1).*sqrt(C.T/C.Rho(1)).*P0(1,:);

P = P0;
V = V0;

%% Simulate interaction
for Iternum = 1:C.Itermax
k1 = dPdT(X, P);
k2 = dPdT(X, P+C.DT.*k1./2);
k3 = dPdT(X, P+C.DT.*k2./2);
k4 = dPdT(X, P+C.DT.*k3);
V = V + (k1+2.*k2+2.*k3+k4)./6;
P = P + V.*C.DT;

%% Plot
if mod(Iternum, 100)==0
PlotStrings(X, P, V);
%        saveas(h, sprintf('Figure%010.f.png', Iternum))
end
end

%% Plot function
function PlotStrings(X, P, V)
figure(1); clf;
C = Constants();

hold on
for StringNum=1:C.N
plot(X(StringNum, :)+StringNum.*max(X, [], 'all'), P(StringNum, :))
yline(0, '--')
end
hold off
ylim([-1 ,1])
end

function output = dPdT(X, P)
C = Constants();

ddPddT = zeros(size(X));

% acceleration
DX = diff(X, 1, 2); DX = DX(:, 1:end-1);
ddPddX = P(:, 1:end-2) - 2.*P(:, 2:end-1) + P(:, 3:end);
ddPddX = ddPddX./DX.^2;
for StringNum=1:C.N
ddPddT(StringNum, 2:end-1) = C.T./C.Rho(StringNum).*ddPddX(StringNum, :);
end

% connections
ddPddT(1:end-1, end) = -C.k.*(P(1:end-1, end) - P(2:end, 1)) +...
C.T.*(P(1:end-1, end-1) - P(1:end-1, end))./(X(1:end-1, end) - X(1:end-1, end-1));
ddPddT(2:end  , 1  ) = -C.k.*(P(2:end, 1) - P(1:end-1, end)) +...
C.T.*(P(2:end  , 2  ) - P(2:end  , 1    ))./(X(2:end  , 2  ) - X(2:end  , 1    ));

% velocity
output = ddPddT.*C.DT;
end


This also requires a file name Constants.m with:

classdef Constants
properties( Constant = true )
Itermax = 5e5;
T = 1;
Rho = [1, 1e3];
N = 1;
K = 100;
DT = 1e-3;
k = 0;
end
end

• Try out periodic boundary conditions ($p_{n+1} \equiv p_1, p_0 \equiv p_n$) to see whether the oscillations are due to reflecting boundary conditions or due to the difference equation itself. May 15 '21 at 14:39
• The phenomenon you observe is known as dispersion. Waves with different wavelength travel with different speeds in the numerical scheme. This properties of numerical methods are studied with dispersion analysis May 15 '21 at 16:35
• Well, numerical dispersion seems to be my problem. At least, It looks very much like this. Do you maybe have a recommended way to solve it? I found some material online, but I would prefer to have a specific direction to focus on. Thanks! May 16 '21 at 22:35

You want to solve $$\ddot P=F(P)$$. In the RK4 stages, you treat the updates as if you were to solve the first order DE $$\dot P=F(P)$$, using the $$k_i$$ values as updates for $$P$$. However in the conclusion you then use the $$k_i$$ as updates for $$V$$ and use a first order step for $$P$$. There is no reason why that should give a recognizable result in any way. At best this all is first order correct. But only if you multiply the derivative with $$DT$$ only once.