# Confusion about Zabusky and Kruskal's stepper for the KdV equation

In Zabusky and Kruskal's paper about solitons, they derive the following update for the Korteweg de Vries equation (their footnote 6): \begin{align*} u_{i}^{j+1} = u_{i}^{j-1} - \frac{1}{3} \frac{k}{h} ( u_{i-1}^{j} + u_{i}^{j} + u_{i+1}^{j})(u_{i+1}^{j} - u_{i-1}^{j}) - \frac{\delta^2k}{h^3} (u_{i+2}^{j} -2u_{i+1}^{j} +2 u_{i-1}^{j} - u_{i-2}^{j}) \end{align*} Their setup is an initial value problem (so that $$u_{i}^{0}$$ is known for all $$i$$) subject to periodic boundary conditions in $$i$$, so that $$u_{i+2N}^{j} = u_{i}^{j}$$. So I have no problem determining the $$i$$ indexes.

But if I only know $$u_{i}^{0}$$, then the update requires $$u_{i}^{1}$$ to get $$u_{i}^{2}$$. (Of course, once I know these, I can get going without problem.) So how can I use this scheme to get $$u_{i}^{1}$$?

(My current suspicion is that they used equation (6) $$u = \cos(\pi(x-ut))$$ for $$u^1$$. Still not sure.)

• Didn't they use a ghost node for $u^{j - 1}$? – nicoguaro Sep 13 '19 at 23:22
• Another option is to use a high-accuracy single step explicit scheme to get $u^1$, then proceed using the multistep iteration. This is what is used in modern multistep integrators. Doing a single time step of RK4, for instance, would probably be more than accurate enough for your needs, and it also conserves momentum – whpowell96 Sep 14 '19 at 3:07
• @nicoguaro: I didn't see anything about a ghost node in their paper, though for my purpose it doesn't matter if I strictly follow exactly their procedure. You will have to forgive my ignorance, but I'm not familiar with how to populate the ghost nodes. – user14717 Sep 14 '19 at 13:36

## 1 Answer

$$\newcommand{\uu}{\mathbf{u}}$$ The method can be condensed to $$\uu^{j+1}=\uu^{j-1}+2kF(\uu^j)$$ which is the weakly stable 2-step Nyström scheme, also called "central Euler method". See some remarks on the scalar method to see in the plot of the error profiles at the end that the method is rather sensitive on the choice of $$\uu^1$$. While the method is second order, an initialization of at least third order and sufficiently small time step is advisable. For simplicity, compute an RK4 step for $$\dot \uu=F(\uu)$$.

However, using the Nyström method or the lsoda solver, both methods produce a mountain scape that corresponds well with the Figure 2 in the paper.

$$N=200$$ sub-intervals in x-diretion, i.e., $$h=0.01$$. Using $$k=0.05\frac{h^3}{\delta^2}$$, with $$\delta=0.022$$from the paper