Numerical instability in finite difference string simulation

Question

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
    % advance time (RK4)
    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

Answer 1

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.

See https://scicomp.stackexchange.com/questions/34257/solving-coupled-odes-using-runge-kutta-method for how the RK4 can be implemented for second order DE, and https://stackoverflow.com/questions/60405185/is-there-a-better-way for what a compact implementation could look like, https://stackoverflow.com/questions/33741338/right-runge-kutta-4 again on the second order to first order system step.