I'm studying a fairly typical problem: a chain of $n$ coupled, non-linear oscillators. Since I want to look at open boundary conditions, the equations of motion for the position of the first and the last oscillator are specified separately:
\begin{align} \ddot{x}_1(t) &= -(x_1(t) - x_2(t)) - V(x_1(t)) + f(t) \\ \ddot{x}_n(t) &= -(x_n(t) - x_{n-1}(t)) - V(x_n(t)) \end{align} where $V(x(t))$ specifies the on-site non-linearity and $f(t)$ is an external driving term. The remaining equations of motion are:
\begin{equation} \ddot{x}_i(t) = -(2 x_i(t) - x_{i+1}(t) - x_{i-1}(t)) - V(x_i(t)), \quad i=2,\dots,n-1 \end{equation}
This is the simplest version of a more general problem I'm trying to understand, but I want to first see how to most efficiently simulate this problem numerically using Mathematica. I've seen many papers where such problems are solved using Molecular Dynamics (MD) simulations where the equations of motions are solved using a Verlet integration algorithm. See for example Sec. II B of https://arxiv.org/abs/0704.1453
Similar to that paper, I want to be able to solve these equations numerically for $n \sim 500$ and for a simulation time $T \sim 1000-5000$, but I'm not sure whether the optimal way to proceed is by using NDSolve or by writing a Verlet algorithm. The code for both methods follows below:
Method 1: Verlet Integration
Brief description of the Verlet Algorithm: a second order differential equation $$\ddot x(t) = F(x(t))$$ with initial conditions $x(0) = X_0$ and $x'(0) = v_0$, can be discretized and numerically solved by this algorithm. First, we pick a time-step $\Delta t$ and define $x_n = x(t_n = n \Delta t)$. Then, the second derivative is approximated as $$ \frac{\Delta^2 x_n}{\Delta t^2} = \frac{x_{n+1} - 2 x_n + x_{n-1}}{\Delta t^2} $$ so that $$ x_{n+1} = 2 x_n - x_{n-1} + \Delta t^2 F(x_n). $$ So in order to find the solution by numerical integration, we set $x_0 = X_0$, $x_1 = X_0 + v_0 \Delta t + \frac{1}{2} \Delta t^2 F(x_0)$, and then iterate $$ x_{i+1} = 2 x_i - x_{i-1} + \Delta t^2 F(x_o), \quad i=1,\dots,n-1. $$
(*Intialize Parameters*)
n = 50; (*Number of Oscillators*)
Tmin = 0; (*Start time*)
Tsim = 100; (*End time*)
tstep = 2000; (*Number of iterations/time-steps*)
h = N[(Tsim - Tmin)/tstep]; (*Time step*)
V[r_] = r^3; (*On-site potential *)
F = 10; (*Drive amplitude*)
\[Omega] = 2.5; (*Drive frequency*)
f[t_] = F Cos[\[Omega] t]; (*Driving term*)
(*Specify Initial Conditions*)
X0 = 0; (*Initial Position*)
V0 = 0; (*Initial Velocity*)
(*Verlet Integration*)
Do[X[i][1] = X0, {i, 1, n}]; (*Set initial positions*)
X[1][2] = X0 + h V0 + h^2/2 F; (*Second step for first oscillator*)
Do[X[i][2] = X0 + h V0 , {i, 2, n}]; (*Second step for remaining oscillators*)
Do[{
X[1][j + 1] = 2 X[1][j] - X[1][j - 1] - h^2 (X[1][j] - X[2][j] - f[(j-1)h] + V[X[1][j]]), (*First Oscillator*)
X[n][j + 1] = 2 X[n][j] - X[n][j - 1] - h^2 (X[n][j] - X[n - 1][j] + V[X[n][j]]), (*Last Oscillator*)
X[i][j + 1] = 2 X[i][j] - X[i][j - 1] - h^2 (2 X[i][j] - X[i - 1][j] - X[i + 1][j] + V[X[i][j]]) (*Remaining Oscillators*)
}, {j, 2, tstep}, {i, 2, n - 1}];
(*Store position data*)
Do[Xdata[i] = Join[{X[i][1], X[i][2]}, Table[X[i][j], {j, 3, tstep + 1}]],{i, 1, n}];
tdata = Table[t, {t, Tmin, Tsim, h}];
Do[Posdata[i] = Transpose[{tdata, Xdata[i]}], {i, 1, n}];
(*Plot Position for i^th oscillator*)
PlotPos[i_] := ListLinePlot[Posdata[i], AxesLabel -> {"t", "y"}, PlotRange -> All]
Method 2: Using NDSolve
(*Intialize Parameters*)
n = 50; (*Number of Oscillators*)
Tmin = 0; (*Start time*)
Tsim = 100; (*End time*)
V[r_] = r^3;(*On-site potential *)
F = 20; (*Drive amplitude*)
\[Omega] = 6; (*Drive frequency*)
f[t_] = F Cos[\[Omega] t]; (*Driving term*)
(*Specify Initial Conditions*)
X0 = 0; (*Initial Position*)
V0 = 0; (*Initial Velocity*)
XN[t_] = Table[Symbol["x" <> ToString[i]][t], {i, 1, n}];
(*Equations of Motion*)
EoM[1] := XN''[t][[1]] - f[t] + (XN[t][[1]] - XN[t][[2]]) + V[XN[t][[1]]] (*First Oscillator*)
EoM[n] := XN''[t][[n]] + (XN[t][[n]] - XN[t][[n - 1]]) + V[XN[t][[n]]](*Last Oscillator*)
EoM[i_] := XN''[t][[i]] + (XN[t][[i]] - XN[t][[i - 1]]) + (XN[t][[i]] - XN[t][[i + 1]]) + V[XN[t][[i]]] (*Remaining Oscillators*)
sol = NDSolve[ArrayFlatten[{Table[EoM[i] == 0, {i, 1, n}], Table[XN[0][[i]] == 0, {i, 1, n}], Table[XN'[0][[i]] == 0, {i, 1, n}]}, 1], XN[t], {t, Tmin, Tsim}];
Comparison
As shown in this plot, both methods give the same solutions:
The first method takes $3.89761$ while the second runs in just $0.17595$ for the same parameters. Using NDSolve is clearly much faster, so I'm wondering if it's better to stick with that or if the MD simulation can be improved to be more efficient, as my algorithm is far from optimized. Even for $n=50$ and $T = 100$, which is much smaller than the parameters I'd like to reach, the Verlet algorithm is taking a long time.
It seems it can be made much better, as in this earlier post: Simulating molecular dynamics efficiently so it would be great if similar speed-up could be achieved for my problem. And if in-built methods are better, then I'm confused as to why people use MD simulations for such problems?
Using either NDSolve or MD simulations, I would appreciate input on how best to proceed to solve this set of equations numerically for large numbers of oscillators and for large simulation times.
UPDATE:
I followed Michael and Henrik's advice for improving my solver by compiling everything. On my system (only 2 cores), my updated code works $\sim 7$ times faster than Alex's optimized Verlet algorithm. See below for comparison:
Alex's Optimized Verlet (I modified the $M$ matrix slightly for open boundary conditions)
n = 64; tmax = 1000; \[Epsilon] = 1.0; m = 1.0; \[Lambda] = \1.0;
x0 = Table[0., {n}]; v0 = Table[0., {n}];
V[x_] := m x + \[Lambda] x^3;
M = SparseArray[{{1, 1} -> -\[Epsilon], {n, n} -> -\[Epsilon], Band[{1, 1}]-> - 2 \[Epsilon], Band[{2, 1}] -> \[Epsilon], Band[{1, 2}] -> \[Epsilon]}, {n, n}]; (*Matrix of Interactions*)
x[t_] = Table[Symbol["x" <> ToString[i]][t], {i, 1, n}];
force[t_] := Table[If[i == 1, 10 Cos[5 t/2], 0], {i, 1, n}];
xN[t_] := Table[Symbol["xN" <> ToString[i]][t], {i, 1, n}];
v[t_] := Table[Symbol["v" <> ToString[i]][t], {i, 1, n}];
h = 1/10; xN[0] = x0; v[0] = v0;
ParallelDo[t1 = t + h;
xN[t1] = xN[t] + v[t] h + (M.xN[t] + force[t] - V[xN[t]]) h^2/2;
v[t1] = v[t] + (M.xN[t] + force[t] - V[xN[t]] + M.xN[t1] + force[t1] -V[xN[t1]]) h/2;, {t, 0, tmax - h, h}]; // AbsoluteTiming
(*{1.39245, Null}*)
My updated algorithm:
VerletDrive[fn_] := With[{f = fn},
Compile[{{x0, _Real, 1}, {v0, _Real, 1}, {tf, _Real}, {n, _Integer}},
Module[{dt, state, Nsim},
Nsim = Length@x0;
dt = N[tf/n];
pos = ConstantArray[0., {n, Nsim}];
pos[[1]] = x0;
pos[[2]] = x0 + dt v0 + 0.5*dt^2*F*ForceArray;
Do[pos[[nc + 1]] = 2 pos[[nc]] - pos[[nc - 1]] + (dt^2)*f@pos[[nc]] + (dt^2)*F*ForceArray*Cos[\[Omega] dt (nc - 1)], {nc, 2, n - 1}];
pos]]]
Nsim = 64; \[Epsilon] = 1.0; m = 1.0; \[Lambda] = 1.0; F = 10.0; \[Omega] = 2.5;
ForceArray = UnitVector[Nsim, 1];
f[x_] := M.x - m x - \[Lambda] x^3
runVer = VerletDrive[f];
x0 = ConstantArray[0.0, Nsim];
v0 = ConstantArray[0.0, Nsim];
tmax = 1000;
h = 0.1;
AbsoluteTiming[data1 = runVer[x0, v0, 10^3, N[tmax/h]];]
(*{0.246074, Null}*)
Comparison during the last 100 seconds:
Compileeverything. $\endgroup$