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

Question

Ι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:

$$U_n(0) = U_{n,0} = a \text{sech}\left(-\frac{L}{2} + nh\right) \ (1) $$

where $ h = \frac{L}{K+1} $ is the distance between the nodes, and $ a > 0 $ is the amplitude. We assume zero initial velocity

$$\dot{U}_n(0) = U_{n,1} = 0.$$

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.

Figure 2.13 shows the dynamics of the initial condition (1) for parameters: $ \omega_0^2 = 1 $, $ K = 99 $, $ L = 200 $, and $\beta = 1 $. 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 $ \delta = 0.1, 0.05, 0.01 $. For $ \delta = 0.1 $, the convergence takes place at the equilibrium point $ \Phi_1 \in C_1 $ which is the basic state.

I want to create the red and the blue plots

When I try to create the first plot (red) I receive this which is correct

(*Define parameters*)\[Omega]d = 1;
k = 99; 
L = 200;
a = 1; 
h = L/(k + 1);
grid = Range[-L/2, L/2, h]; (*Spatial points*)
kappa = Length[grid]; (*number of spatial points*)

(*Initial conditions*)
initialDisplacement[n_] := a Sech[-L/2 + n h];
initialVelocity[n_] := 0;

(*Prepare the data points*)
dataPoints = Table[{grid[[n]], initialDisplacement[n]}, {n, 1, kappa}];

(*Plot the initial condition*)
initialPlot = 
  ListPlot[dataPoints, PlotRange -> {{-L/2, L/2}, {-2, 2}}, 
   AxesLabel -> {"x_n", "U_n"}, PlotLabel -> "t=0", Joined -> True, 
   PlotStyle -> {Red}, 
   Epilog -> {Red, PointSize[Medium], Point[dataPoints]}, 
   FrameLabel -> {"\!\(\*TemplateBox[<|\"boxes\" -> FormBox[\n\
SubscriptBox[\nStyleBox[\"x\", \"TI\"], \nStyleBox[\"n\", \"TI\"]], \
TraditionalForm], \"errors\" -> {}, \"input\" -> \"x_n\", \"state\" -> \
\"Boxes\"|>,\n\"TeXAssistantTemplate\"]\)", 
     "\!\(\*TemplateBox[<|\"boxes\" -> FormBox[\nSubscriptBox[\n\
StyleBox[\"u\", \"TI\"], \nStyleBox[\"n\", \"TI\"]], \
TraditionalForm], \"errors\" -> {}, \"input\" -> \"u_n\", \"state\" -> \
\"Boxes\"|>,\n\"TeXAssistantTemplate\"]\)"}, 
   PlotTheme -> "Scientific", 
   Ticks -> {{-100, -50, 0, 50, 100}, {-2, -1, 0, 1, 2}}, 
   AxesOrigin -> {0, 0}];

(*Display the plot*)
initialPlot

Then I tried the following for different values of $\delta$. What should I change to my code?

(* Parameters *)
a = 2;
deltaList = {0.1, 0.05, 0.01};
L = 200;
omegad = 1;
beta = 1;
h = 2;
grid = Range[-L/2, L/2, h];
k = Length[grid];
omegad2 = (h^2) * omegad^2;
tmax = 3000;

(* Function to solve and plot the system for a given delta *)
solveAndPlot[delta_] := Module[{var, eqs, bc, ic, sol, initialPlot, finalPlot, deltaHat},
  deltaHat = h * delta;
  var = Table[u[i][t], {i, k}];
  
  eqs = Table[
    u[i]''[t] - (u[i + 1][t] + u[i - 1][t] - 2 u[i][t]) + 
    deltaHat u[i]'[t] - omegad2 (u[i][t] - beta u[i][t]^3) == 0, 
    {i, 2, k - 1}
  ];
  
  bc = {u[1][t] == 0, u[k][t] == 0};
  
  ic = Join[
    Thread[var == a Sech[-L/2 + grid]] /. t -> 0, 
    Thread[D[var, t] == 0] /. t -> 0
  ];
  
  sol = NDSolve[Join[eqs, bc, ic], var, {t, 0, tmax}];
  
  initialPlot = ListPlot[
    Transpose[{grid, var /. sol[[1]] /. t -> 0}], 
    FrameLabel -> {"x_n", "u_n"}, 
    PlotLabel -> Style["t = 0", Italic], 
    PlotTheme -> "Scientific", PlotStyle -> {Red, PointSize[Medium]}
  ];
  
  finalPlot = ListPlot[
    Transpose[{grid, var /. sol[[1]] /. t -> tmax}], 
    FrameLabel -> {"x_n", "u_n"}, 
    PlotLabel -> Style["t = 3000, delta = " <> ToString[delta], Italic], 
    PlotTheme -> "Scientific", PlotStyle -> {Blue, PointSize[Medium]}
  ];
  
  {initialPlot, finalPlot}
];

(* Generate and display plots for each delta *)
plots = Flatten[Table[solveAndPlot[delta], {delta, deltaList}], 1];
GraphicsGrid[Partition[plots, 2], ImageSize -> Large]

Answer 1

Initial condition should be symmetric around 0, while in your code there is some shift. For symmetric condition we have

(*parameters*)a = 1;(*Amplitude of the initial condition*)
j = 2; (*Mode number*)
L = 200; (*Length of the system*)

omegad = 1; (*characteristic frequency*)
beta = 1;(*Nonlinearity parameter*)

(*damping coefficient*)
delta = 0.05;

(*spatial grid*)
h = 1;
grid = Range[-L/2, L/2, h];(*Spatial points*)
k = Length[grid]; (*number of spatial points*)
(*discretization parameter*)
(*Scaled parameters*)
omegad2 = (h^2)*omegad^2;
deltaHat = h*delta;

tmax = 3000; (*Maximum time*)


(*Initial positions*)

var = Table[u[i][t], {i, k}]; eqs = 
 Table[u[i]''[t] - (u[i + 1][t] + u[i - 1][t] - 2 u[i][t]) + 
    deltaHat u[i]'[t] - omegad2 (u[i][t] - beta u[i][t]^3) == 0, {i, 
   2, k - 1}]; bc = {u[1][t] == 0, u[k][t] == 0};
ic = Join[Thread[var == a  Sech[grid]] /. t -> 0, 
   Table[u[i]'[0] == 0, {i, 2, k - 1}]];
sol = NDSolve[Join[eqs, {u[1]'[t] == 0, u[k]'[t] == 0}, ic], 
    var, {t, 0, tmax}];

Visualization

{Show[ListPlot[Transpose[{grid, var /. sol[[1]] /. t -> 0}], 
   FrameLabel -> {"\!\(\*SubscriptBox[\(X\), \(n\)]\)", 
     "\!\(\*SubscriptBox[\(U\), \(n\)]\)"}, 
   PlotLabel -> "Initial Positions of DNA Segment at t = 0", 
   PlotTheme -> "Scientific", PlotStyle -> {Red, PointSize[Medium]}, 
   PlotRange -> All], 
  ListLinePlot[Transpose[{grid, var /. sol[[1]] /. t -> 0}], 
   PlotStyle -> {Red, Dashed}, PlotRange -> All]], 
 Show[ListPlot[Transpose[{grid, var /. sol[[1]] /. t -> tmax}], 
   FrameLabel -> {"\!\(\*SubscriptBox[\(X\), \(n\)]\)", 
     "\!\(\*SubscriptBox[\(U\), \(n\)]\)"}, 
   PlotLabel -> "Positions of DNA Segment at t = 3000", 
   PlotTheme -> "Scientific", PlotStyle -> {Blue, PointSize[Medium]}, 
   PlotRange -> All], 
  ListLinePlot[Transpose[{grid, var /. sol[[1]] /. t -> tmax}], 
   PlotStyle -> {Blue, Dashed}]]}