Discretisation

The simplest discretisation of (2.61), based on the Euler method, gives the equation

$$y_j^{(n+1)} = y_j^{(n)} - \left[{\mathchoice{{\textstyle{\f... ...1}^{(n)} + 2 y_{j-1}^{(n)} - y_{j-2}^{(n)}\right\}\right].$$ (2.63)

You should check that this is indeed a sensible discretisation of (2.61) but that it is unstable. Note that when analysing the non-linear term in (2.63) you should make the substitution $y\mapsto y + \delta y$ and retain the linear terms in $\delta y$. Thereafter you should treat $y$ as a constant, independent of $x$ and $t$, and apply the von Neumann method to $\delta y$. If you find the full stability analysis difficult you might consider the 2 limits of large and small $y$. In the former case the 3rd derivative is negligible and (2.61) reduces to a non-linear advection equation, whereas in the latter the non-linear term is negligible and the equation is similar to the diffusion equation but with a 3rd derivative. In any case you will require to choose $\delta t$ so that the equation is stable in both limits. You are free to choose any method you wish to solve the equation but you will find the Runge-Kutta or Predictor-Corrector methods most reliable. Hence treating $y_j$ as a long vector $\bi{y}$ and the terms on the right-hand-side of (2.63) as a vector function $\mathop{\bf f}\nolimits(\bi{y})$ the R-K method can be written concisely as

$$\bi{y}'^{(n+{\mathchoice{{\textstyle{\frac12}}}{{\textstyle{\frac12}}}{{\scriptstyle{1/2}}}{{\scriptscriptstyle{1/2}}}})} = \bi{y}^{(n)} - {\mathchoice{{\textstyle{\frac12}}}{{\textstyle{\f... .../2}}}{{\scriptscriptstyle{1/2}}}}\delta t \mathop{\bf f}\nolimits(\bi{y}^{(n)})$$ (2.64)

$$\bi{y}^{(n+1)} = \bi{y}^{(n)} - \delta t \mathop{\bf f}\nolimits(\bi{y}'^{(n+{\mat... ...\textstyle{\frac12}}}{{\scriptstyle{1/2}}}{{\scriptscriptstyle{1/2}}}})})\quad,$$ (2.65)

where $\delta t\mathop{\bf f}\nolimits(\bi{y})$ is the quantity in square brackets $[]$ in (2.63). Check that this method is stable, at least in the 2 limiting cases. Bear in mind that the Runge-Kutta method is usable for oscillatory equations in spite of the small instability as long as the term $(\omega\delta t)^4$ is small. By studying the analytical solution (2.62) you should be able to choose a sensible value for $\delta x$ in terms of $\alpha$ and from the stability conditions you can deduce an appropriate $\delta t$. Again, by looking at (2.62) you should be able to decide on a sensible size for the total system. You should use periodic boundary conditions, so that your solitons can run around your system several times if necessary. The easiest way to do this is by using ``ghost'' elements at each end of your arrays. Suppose your arrays should run from $y_1$ to $y_N$. Then you can add a couple of extra elements to each end: $y_{-1}, y_0, y_{N+1}, y_{N+2}$. After each step you can then assign these values as

$$y_{-1} = y_{N-1},\quad y_0 = y_N,\quad y_{N+1} = y_1, \quad y_{N+2} = y_2$$ (2.66)

so that the derivative terms in (2.63) can be calculated without having to take any special measures.