Numerical Simulation of a Damped, Driven Nonlinear Wave System with Spatially Extended Initial Conditions

The study of the dynamics of the discrete Klein - Gordon equation (DKG) with friction is given by the equation :

In the above equation, W describes the potential function:

to which every coupled unit adheres. In Eq. (1), the variable is the unknown displacement of the oscillator occupying the n-th position of the lattice, and is the discretization parameter. We denote by h the distance between the oscillators of the lattice. The chain (DKG) contains linear damping with a damping coefficient , while β is the coefficient of the nonlinear cubic term.

For the DKG chain (1), we will consider the problem of initial-boundary values, with initial conditions

and Dirichlet boundary conditions at the boundary points and , that is,

.

Therefore, when necessary, we will use the short notation for the one-dimensional discrete Laplacian

We investigate numerically the dynamics of the system (1)-(2)-(3). Our first aim is to conduct a numerical study of the property of Dynamic Stability of the system, which directly depends on the existence and linear stability of the branches of equilibrium points.

For the discussion of numerical results, it is also important to emphasize the role of the parameter . By changing the time variable , we rewrite Eq. (1) in the form

The change in the scaling of the lattice parameter of the problem makes it important to comment here on the nature of the continuous and anti-continuous limit. For the anti-continuous limit, we need to consider the case in Eq. (1) where . In the case of nonlinearity present in the governing equations, the continuous limit needs to be taken as well. On the other hand, for small values of the parameter, the system becomes significant. However, because of the model we consider, we take the asymptotic linear limit as .

We consider spatially extended initial conditions of the form:

where is the distance of the grid and is the amplitude of the initial condition

We also assume zero initial velocity:

I am trying to create the following plots but as you can see my code doesn;t give these results. Any sugestions?

• for

% Parameters
a = 2;         % Amplitude of the initial condition
j = 2;         % Mode number
L = 200;       % Length of the system
K = 99;       % Number of spatial points

% Spatial grid
h = L / (K + 1);
n = -L/2:h:L/2;  % Spatial points


% Compute U_n(0) for each n
U_n0 = a * sin((j * pi * h * n) / L);

% Plot the results
figure;
plot(n, U_n0, 'ro'); % 'ro' creates red circles
xlabel('$x_n$', 'Interpreter', 'latex');
ylabel('$U_n$', 'Interpreter', 'latex');
title('$t=0$', 'Interpreter', 'latex');
xlim([-L/2 L/2]);
ylim([-3 3]);
grid on;

UPDATED : I need to replicate the red and the blue graphs of the papers. I am sure that I need to change a parameter but I dont know which one for and .

% Parameters
L = 200;  % Length of the system
K = 99;  % Number of spatial points
j = 2;  % Mode number
omega_d = 1;  % Characteristic frequency
beta = 1;  % Nonlinearity parameter
delta = 0.05;  % Damping coefficient

% Spatial grid
h = L / (K + 1);
n = linspace(-L/2, L/2, K+2);  % Spatial points
N = length(n);
omegaDScaled = h * omega_d;
deltaScaled = h * delta;

% Time parameters
dt = 0.05; % Time step
tmax = 3000; % Maximum time
tspan = 0:dt:tmax; % Time vector

% Values of amplitude 'a' to iterate over
a_values = [2, 1.95, 1.9, 1.85, 1.82];  % Modify this array as needed

% Differential equation solver function
function dYdt = odefun(~, Y, N, h, omegaDScaled, deltaScaled, beta)
    U = Y(1:N);
    Udot = Y(N+1:end);
    Uddot = zeros(size(U));

    % Laplacian (discrete second derivative)
    for k = 2:N-1
        Uddot(k) = (U(k+1) - 2 * U(k) + U(k-1)) / h^2;
    end

    % System of equations
    dUdt = Udot;
    dUdotdt = Uddot - deltaScaled * Udot + omegaDScaled^2 * (U - beta * U.^3);

    % Pack derivatives
    dYdt = [dUdt; dUdotdt];
end

% Create a figure for subplots
figure;

% Initial plot
a_init = 2;  % Example initial amplitude for the initial condition plot
U0_init = a_init *  sin((j * pi * h * n) / L); % Initial displacement
U0_init(1) = 0; % Boundary condition at n = 0
U0_init(end) = 0; % Boundary condition at n = K+1
subplot(3, 2, 1);
plot(n, U0_init, 'r.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
xlabel('$x_n$', 'Interpreter', 'latex');
ylabel('$U_n$', 'Interpreter', 'latex');
title('$t=0$', 'Interpreter', 'latex');
set(gca, 'FontSize', 12, 'FontName', 'Times');
 xlim([-L/2 L/2]);
ylim([-3 3]);
grid on;

% Loop through each value of 'a' and generate the plot
for i = 1:length(a_values)
    a = a_values(i);

    % Initial conditions
    U0 = a * sin((j * pi * h * n) / L); % Initial displacement
    U0(1) = 0; % Boundary condition at n = 0
    U0(end) = 0; % Boundary condition at n = K+1
    Udot0 = zeros(size(U0)); % Initial velocity

    % Pack initial conditions
    Y0 = [U0, Udot0];

    % Solve ODE
    opts = odeset('RelTol', 1e-5, 'AbsTol', 1e-6);
    [t, Y] = ode45(@(t, Y) odefun(t, Y, N, h, omegaDScaled, deltaScaled, beta), tspan, Y0, opts);

    % Extract solutions
    U = Y(:, 1:N);
    Udot = Y(:, N+1:end);

    % Plot final displacement profile
    subplot(3, 2, i+1);
    plot(n, U(end,:), 'b.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
    xlabel('$x_n$', 'Interpreter', 'latex');
    ylabel('$U_n$', 'Interpreter', 'latex');
    title(['$t=3000$, $a=', num2str(a), '$'], 'Interpreter', 'latex');
    set(gca, 'FontSize', 12, 'FontName', 'Times');
      xlim([-L/2 L/2]);
ylim([-2 2]);
    grid on;
end

% Adjust layout
set(gcf, 'Position', [100, 100, 1200, 900]); % Adjust figure size as needed

[William Rose]

@Athanasios Paraskevopoulos,

In your first equation, does

represent

?

[Athanasios Paraskevopoulos]

Exactly

[William Rose]

@Athanasios Paraskevopoulos, I'm sorry, but I am still not understanding things. I assume "n" represents spatial location. What is the difference between indices j and n, in the equation for the initial condition for position:

It seems to me that the initial condition should be

where n=1...K. I don;t understand the notation. Also, you write

which makes me think the subscript 0 in indicates time=0. But the initital condition for velocity is

where, according to my earlier interpretation, would represent the position at location=n and time t=1. I guess I don't understand the notation.

[William Rose]

@Athanasios Paraskevopoulos, I think my equation above for is incorrect, because I forgot to divide by h^2. Shouldn't it be

?

[William Rose]

@Athanasios Paraskevopoulos,

[edit: fix spelling errors]

I also want to be sure I understand the units. I assume U is position. Let's use meters for position. Then the units are and . What are the units for δ? Whatever those units are will determine the units for . (I assume h has units of length.)

The left hand side of your first equation is

The units for the first two terms are incompatible. If has units of length, which seems quite possible, then the third term units do not equal the units for the first or second term.

The right hand side is:

Therefore β must have units of 1/length = 1/meters. Does that seem right to you? What is the physical meanng of beta in this problem? I assume L is a length, and therefore h is a length. Then has units length^2/time^2. But if that is true, then the units for the right hand side are m^3/s^2, which is incompatible with the left hand side.

Please help me understand the units.

[William Rose]

@Athanasios Paraskevopoulos, it looks like beads on a string. The simple linear equation for beads on a string is

where (without h^2 in the denominator), m=bead mass, h=distance between beads, T=tension in the string (force). This equation assumes the displacements are small: therefore the magnitude (but not direction) of tension is constant at all times and places, and sin(theta)~=theta. The constant T/(mh) makes the units work out correctly.

If the beads are moving in a viscous fluid, with friction proportional to velocity, the equation of motion is

where , where η=viscosity (N-s/m) and m=particle mass. You define , where h=distance between beads.

Your equation's left side equals what is shown above, except for the constant T/(mh). On the right side you have . I am not sure what the physical meaning of this is. I think it would be unstable if beta=0.

[Athanasios Paraskevopoulos]

The given equation and the initial conditions describe the dynamics of a spatially discrete system, often encountered in the context of lattice vibrations or discrete models in physics. Let's break down the notation and terms step by step.

• The Equation

The primary equation is:

• Terms in the Equation:

- : The displacement at site and time t.

- : The second time derivative of , representing the acceleration at site n.

- : The first time derivative of, representing the velocity at site n.

- : The discrete Laplacian operator applied to , which typically represents the interaction between neighboring sites. For a one-dimensional lattice, it can be defined as:

- : A damping term proportional to the velocity, with being the damping coefficient.

- : A squared frequency term, defined as , where is a spatial scaling factor, and is a characteristic frequency.

- A nonlinear coefficient that dictates the strength of the nonlinear term

• Parameters

- : The squared frequency term incorporating the spatial scaling factor h.

- : The damping coefficient incorporating the spatial scaling factor

• Initial Conditions

The initial conditions are given for and its time derivative:

1. Initial Displacement:

- a: The amplitude of the initial displacement.

- : A sinusoidal function representing the initial displacement pattern across the lattice.

- : The grid spacing, where L is the length of the domain andK is the number of segments.

2. Initial Velocity:

- This indicates that the initial velocity of all sites is zero.

[William Rose]

@Athanasios Paraskevopoulos, thank you. More later.

[Athanasios Paraskevopoulos]

@William Rose I tried to explain it better in my initial post

[William Rose]

@Athanasios Paraskevopoulos, thanhk you for your expanded explanation.

Your revised initial post begins with

and in the following text you say " is the discretization parameter". Are K and k the same thing?

The left had side of equation 1 includes W'(Un), but you define , without the prime after W. What is the meaning of the prime after W in equation 1? Are W'() and W() the same thing?

In an earlier comment, I asked about consistency of units. I ask about units because it is a good idea to check for consistent units, when a simulation is not working as expected - which is what you have reported. Your revised post does not address those questions, so I wil ask th questions again, in a revised form, in light of your revised initial post:

1. You define a new time variable, t=t/h^2. In my experience, normalized variables are frequenctly dimensionless. For example, time in seconds becomes t/T, which is dimensionless, or x in meters becomes x/L, which is dimensionless. But your new time variable does not work this way. In fact, the new time variable has units of time/length^2. Why make the units more complicated?
2. What are the units for K, in equation 1? I think K=c^2, where c=speed of wave propagation, with units length/time, and therefore units for K are length^2/time^2. Do you agree?
3. What are the units for δ in equation 1? I think the units are 1/time. Do you agree?

I think the answer to my quesiton about new time=t/h^2 is that you are choosing a time variable in which c, the speed of wave propagation (in the absence of damping), is unity. Is that right?

You define

K = 100;  % Number of spatial points
h = L / (K + 1);
n = -L/2:h:L/2;  % Spatial points

This results in values for n that are not integers, and it results in 102 spatial points, where I think you want 100 spatial points. If you use

h = L / (K - 1);

this inconsistency will disappear.

More later.

[Athanasios Paraskevopoulos]

@William Rose Fo the we define as the discretization parameter.

When

So, the derivative is:

I want to create the following graphs for and

[William Rose]

@Athanasios Paraskevopoulos, if W'(U) in equation 1 means d(W(U))/DU, then I would expect the right side of this equation which you posted

to be

and it is not.

[Athanasios Paraskevopoulos]

I am studying the PhD thesis Dynamics of nonlinear lattices: asymptotic behavior and study of the existence and stability of tracked oscillations and I have observed a mistake at a function on page 69. I was trying to prove the result you said. I just wrote the equation as it is in the PhD thesis at this post

[William Rose]

@Athanasios Paraskevopoulos, I will get you somehting on this, but I'm busy for the next two days. Sorry for the delay. And I apologize for asking so many questions.

[Athanasios Paraskevopoulos]

@William Rose Thank you so much! I appreciate it! Don 't worry I am asking also many questions too,

I am feeling that I am getting closer to the solution but also i feel far

Answer 1

[Athanasios Paraskevopoulos]

I found it i should remove the discretization parameter because in the following.

I understood that

@William Rose Thank you so much! Your questions helped me a lot

% Parameters
L = 200;  % Length of the system
K = 99;  % Number of spatial points
j = 2;  % Mode number
omega_d = 1;  % Characteristic frequency
beta = 1;  % Nonlinearity parameter
delta = 0.05;  % Damping coefficient

% Spatial grid
h = L / (K + 1);
n = linspace(-L/2, L/2, K+2);  % Spatial points
N = length(n);
omegaDScaled = h * omega_d;
deltaScaled = h * delta;

% Time parameters
dt = 1; % Time step
tmax = 3000; % Maximum time
tspan = 0:dt:tmax; % Time vector

% Values of amplitude 'a' to iterate over
a_values = [2, 1.95, 1.9, 1.85, 1.82];  % Modify this array as needed

% Differential equation solver function
function dYdt = odefun(~, Y, N, h, omegaDScaled, deltaScaled, beta)
    U = Y(1:N);
    Udot = Y(N+1:end);
    Uddot = zeros(size(U));

    % Laplacian (discrete second derivative)
    for k = 2:N-1
        Uddot(k) = (U(k+1) - 2 * U(k) + U(k-1)) ;
    end

    % System of equations
    dUdt = Udot;
    dUdotdt = Uddot - deltaScaled * Udot + omegaDScaled^2 * (U - beta * U.^3);

    % Pack derivatives
    dYdt = [dUdt; dUdotdt];
end

% Create a figure for subplots
figure;

% Initial plot
a_init = 2;  % Example initial amplitude for the initial condition plot
U0_init = a_init *  sin((j * pi * h * n) / L); % Initial displacement
U0_init(1) = 0; % Boundary condition at n = 0
U0_init(end) = 0; % Boundary condition at n = K+1
subplot(3, 2, 1);
plot(n, U0_init, 'r.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
xlabel('$x_n$', 'Interpreter', 'latex');
ylabel('$U_n$', 'Interpreter', 'latex');
title('$t=0$', 'Interpreter', 'latex');
set(gca, 'FontSize', 12, 'FontName', 'Times');
 xlim([-L/2 L/2]);
ylim([-3 3]);
grid on;

% Loop through each value of 'a' and generate the plot
for i = 1:length(a_values)
    a = a_values(i);

    % Initial conditions
    U0 = a * sin((j * pi * h * n) / L); % Initial displacement
    U0(1) = 0; % Boundary condition at n = 0
    U0(end) = 0; % Boundary condition at n = K+1
    Udot0 = zeros(size(U0)); % Initial velocity

    % Pack initial conditions
    Y0 = [U0, Udot0];

    % Solve ODE
    opts = odeset('RelTol', 1e-5, 'AbsTol', 1e-6);
    [t, Y] = ode45(@(t, Y) odefun(t, Y, N, h, omegaDScaled, deltaScaled, beta), tspan, Y0, opts);

    % Extract solutions
    U = Y(:, 1:N);
    Udot = Y(:, N+1:end);

    % Plot final displacement profile
    subplot(3, 2, i+1);
    plot(n, U(end,:), 'b.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
    xlabel('$x_n$', 'Interpreter', 'latex');
    ylabel('$U_n$', 'Interpreter', 'latex');
    title(['$t=3000$, $a=', num2str(a), '$'], 'Interpreter', 'latex');
    set(gca, 'FontSize', 12, 'FontName', 'Times');
      xlim([-L/2 L/2]);
ylim([-2 2]);
    grid on;
end

% Adjust layout
set(gcf, 'Position', [100, 100, 1200, 900]); % Adjust figure size as needed