Plotting discrete Klein - Gordon equation (DKG) with friction in DNA

I am exploring 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 :

The objective of this simulation is to model the dynamics of a segment of DNA under thermal fluctuations with fixed boundaries using a modified discrete Klein-Gordon equation. The model incorporates elasticity, nonlinearity, and damping to provide insights into the mechanical behavior of DNA under various conditions.

My question to the following code is: Why are the horizontal lines in the first and the last plot? I can not undesrastant that

% Parameters
numBases = 200; % Number of base pairs, representing a segment of DNA
kappa = 0.1; % Elasticity constant
omegaD = 0.2; % Frequency term
beta = 0.05; % Nonlinearity coefficient
delta = 0.01; % Damping coefficient
% Random initial perturbations to simulate thermal fluctuations
initialPositions = 0.01 + (0.02-0.01).*rand(numBases,1);
initialVelocities = zeros(numBases,1); % Assuming initial rest state

% Define the differential equations
dt = 0.05; % Time step
tmax = 50; % Maximum time
tspan = 0:dt:tmax; % Time vector
x = zeros(numBases, length(tspan)); % Displacement matrix
x(:,1) = initialPositions; % Initial positions

% Velocity-Verlet algorithm for numerical integration
for i = 2:length(tspan)
    % Boundary conditions: fixed ends
    x(1, i) = 0;
    x(numBases, i) = 0;
    
    % Compute acceleration
    acceleration = zeros(numBases,1);
    for n = 2:numBases-1
        acceleration(n) = kappa * (x(n+1, i-1) - 2 * x(n, i-1) + x(n-1, i-1)) ...
            - delta * initialVelocities(n) - omegaD^2 * (x(n, i-1) - beta * x(n, i-1)^3);
    end
    
    % Update positions and velocities
    x(:, i) = x(:, i-1) + dt * initialVelocities + 0.5 * dt^2 * acceleration;
    initialVelocities = initialVelocities + dt * acceleration;
end

% Visualization of displacement over time for each base pair
figure;
hold on;
for n = 1:numBases
    plot(tspan, x(n, :));
end
xlabel('Time');
ylabel('Displacement');
legend(arrayfun(@(n) ['Base ' num2str(n)], 1:numBases, 'UniformOutput', false));
title('Displacement of DNA Bases Over Time');
hold off;

1212

% Improved 3D plot for displacement
figure;
[X, T] = meshgrid(1:numBases, tspan);
surf(X', T', x, 'EdgeColor', 'none');
xlabel('Base Pair');
ylabel('Time (s)');
zlabel('Displacement (m)');
title('3D View of DNA Base Displacements');
colormap(jet); % Colorful gradient representing displacement magnitude
shading interp; % Interpolated shading for smooth appearance
colorbar; % Adds a color bar to indicate displacement magnitude

% Snapshot visualization at a specific time
snapshotTime = 10; % Desired time for the snapshot
[~, snapshotIndex] = min(abs(tspan - snapshotTime)); % Find closest index
snapshotSolution = x(:, snapshotIndex); % Extract displacement at the snapshot time

% Plotting the snapshot
figure;
stem(1:numBases, snapshotSolution, 'filled'); % Discrete plot using stem
title(sprintf('DNA Model Displacement at t = %d seconds', snapshotTime));
xlabel('Base Pair Index');
ylabel('Displacement');

% Time vector for detailed sampling
tDetailed = 0:0.5:50; % Detailed time steps

% Initialize an empty array to hold the data
data = [];

% Generate the data for 3D plotting
for i = 1:numBases
    % Interpolate to get detailed solution data for each base pair
    detailedSolution = interp1(tspan, x(i, :), tDetailed);
    % Concatenate the current base pair's data to the main data array
    data = [data; repmat(i, length(tDetailed), 1), tDetailed', detailedSolution'];
end

% 3D Plot
figure;
scatter3(data(:,1), data(:,2), data(:,3), 10, data(:,3), 'filled');
xlabel('Base Pair');
ylabel('Time');
zlabel('Displacement');
title('3D Plot of DNA Base Pair Displacements Over Time');

colorbar; % Adds a color bar to indicate displacement magnitude

Answer 1

[Torsten]

acceleration(1) = acceleration(numBases) = 0 

for every value of i in your loop.

From the equation

initialVelocities = initialVelocities + dt * acceleration;

it follows that

initialVelocities(1) = initialVelocities(numBases) = 0

for every value of i in your loop.

From the equation

x(:, i) = x(:, i-1) + dt * initialVelocities + 0.5 * dt^2 * acceleration;

it follows that x(1,i) and x(numBases,i) remain unchanged and are equal to x(1,1) and x(numBases,1) for every value of i in the loop.

If you want to keep them as 0, you will have to update as

x(2:numBases-1, i) = x(2:numBases-1, i-1) + dt * initialVelocities(2:numBases-1) + 0.5 * dt^2 * acceleration(2:numBases-1);

[Pat Gipper]

Without much insight I tried inserting the following lines of code starting at line 31 as a guess for the 1st and 200th acceleration values.

acceleration(1)=acceleration(2);% Use #2 as a guess
acceleration(200)=acceleration(199);% Use #199 as a guess

[Torsten]

I think position and velocity in the boundary points have to be prescribed.

OP chose to set position = velocity = 0 in both endpoints.

"acceleration" follows from these settings.

[Athanasios Paraskevopoulos]

Thank you that was really helpful