Plotting the Evolution of Spatially Localized Initial Conditions for Discrete Klein-Gordon Equation

Ιn continuation of the study with which I have been fascinated Numerical Simulation of a Damped, Driven Nonlinear Wave System with Spatially Extended Initial Conditions

Now we examine the role of damping in the structure of the branches of equilibrium points. In this study, we consider localized initial conditions of the form:

where is the distance between the nodes, and is the amplitude. We assume zero initial velocity

We also note that if we consider an infinite chain, the initial conditions $(1)$ satisfy the Dirichlet boundary conditions only asymptotically. However, since the smallest half-length of the chain has been set to $ L/2 = 100$, the error is of the order of $ 10^{-44} $, which does not affect numerical computations.

The figure shows the dynamics of the initial condition (1) for parameters: and . In the first image of Figure, we observe the profile of the initial condition.

The subsequent images show the dynamics of the system for different values of the damping force . For , the convergence takes place at the equilibrium point which is the basic state.

I want to create the red and the blue plots

Now I have tried but as you can see the result are different. What should I change?

% Parameters
L = 200;  % Length of the system
K = 99;  % Number of spatial points
omega_d = 1;  % Characteristic frequency
beta = 1;  % Nonlinearity parameter
delta_values = [0.1, 0.05, 0.01];  % Damping coefficients

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

% 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];  % Initial amplitude for the initial condition plot

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

% Loop through each value of 'delta' and generate the plot
for i = 1:length(delta_values)
    delta = delta_values(i);
    deltaScaled = h * delta;
    a = a_values(1);
    % Initial conditions
    U0 = a * sech(-L/2 + n*h); % 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(2, 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$, $\delta=', num2str(delta), '$'], 'Interpreter', 'latex');
    set(gca, 'FontSize', 12, 'FontName', 'Times');
    xlim([-L/2 L/2]);
    ylim([-2 2]);
    grid on;
end

% Initial plot
subplot(2, 2, 1);
a_init = 2;  % Example initial amplitude for the initial condition plot
U0_init = a_init * sech(-L/2 + n*h); % Initial displacement
U0_init(1) = 0; % Boundary condition at n = 0
U0_init(end) = 0; % Boundary condition at n = K+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;

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

Answer 1

[Athanasios Paraskevopoulos]

% Parameters
L = 200;  % Length of the system
K = 99;  % Number of spatial points
omega_d = 1;  % Characteristic frequency
beta = 1;  % Nonlinearity parameter
delta_values = [0.1, 0.05, 0.01];  % Damping coefficients

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

% 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 = [1];  % Initial amplitude for the initial condition plot

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

% Loop through each value of 'delta' and generate the plot
for i = 1:length(delta_values)
    delta = delta_values(i);
    deltaScaled = h * delta;
    a = a_values(1);
    % Initial conditions
    U0 = a * sech( n*h); % 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(2, 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$, $\delta=', num2str(delta), '$'], 'Interpreter', 'latex');
    set(gca, 'FontSize', 12, 'FontName', 'Times');
    xlim([-L/2 L/2]);
    ylim([-2 2]);
    grid on;
end

% Initial plot
subplot(2, 2, 1);
a_init = 1;  % Example initial amplitude for the initial condition plot
U0_init = a_init * sech( n*h); % Initial displacement
U0_init(1) = 0; % Boundary condition at n = 0
U0_init(end) = 0; % Boundary condition at n = K+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;

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