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))

%% 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, '--')
hold off
ylim([-1 ,1])

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, :);

    % 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;

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;
  • $\begingroup$ 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. $\endgroup$
    – uranix
    May 15 '21 at 14:39
  • 1
    $\begingroup$ 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 $\endgroup$
    – uranix
    May 15 '21 at 16:35
  • $\begingroup$ 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! $\endgroup$
    – DeadlosZ
    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.

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.

  • $\begingroup$ Hi, thanks for the heads up! I implemented the RK4 method as illustrated in the link, but the effect still exists and looks about the same. Any ideas? $\endgroup$
    – DeadlosZ
    May 16 '21 at 20:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.