*2.2. Semi-Implicit Finite Difference Method*

The problem (1) defined in (−∞, ∞) can be restricted to a finite interval Ω = [−*b*, *b*] for a positive constant *b*. This Initial Boundary Value Problem has a third order derivative with respect to *x* and, due to the presence of a forcing, the direction of wave propagation is unknown. Therefore, it is not a trivial task to impose the correct boundary conditions at finite boundaries. This leads us to solve the alternative form of the forced KdV equation. This results in the following Initial Boundary Value problem

$$\begin{aligned} u\_{tx}(\mathbf{x},t) + 2\lambda u\_{xx}(\mathbf{x},t) - \frac{3}{2}(\mathbf{u}^2)\_{xx}(\mathbf{x},t) - \frac{1}{3}u\_{xxxx}(\mathbf{x},t) &= f\_{xx}(\mathbf{x}), \\ u(\mathbf{x},0) = u\_0(\mathbf{x}) \\ u\_x(\mathbf{x},t) = u\_{xx}(\mathbf{x},t) &= 0 \text{ at } \mathbf{x} = \pm b, \end{aligned} \tag{3}$$

*Symmetry* **2018**, *10*, 129

The problem involves the fourth derivative for the spatial variable, therefore four boundary conditions are required. This is done by imposing the zero Neumann boundary conditions at each boundary. Now, we employ the finite difference scheme to find the numerical solutions of equation above. Let us denote *x*0 = <sup>−</sup>*b*, *x*1 = −*b* + Δ*x*, ··· , *xi* = −*b* + *i*Δ*x*, ··· , *xN* = *b*, *t*0 = 0, *t*1 = Δ*t*, ··· , *tn* = *<sup>n</sup>*Δ*t*, ··· , *tNt* = *T*, *<sup>u</sup>*(*xi*, *tn*) = *uni* and *fxx*(*xi*) = ¯ *fi* where Δ*x* = 2*b*/*N* and Δ*t* = *T*/*Nt* with positive integers *N* and *Nt*, respectively. First, the Crank-Nicolson scheme with standard difference schemes is applied to (3) and the resulting finite difference equations are written as

$$\frac{D^0(\boldsymbol{u}\_i^{n+1} - \boldsymbol{u}\_i^n)}{2\Delta\mathbf{x}\Delta t} + 2\lambda\Delta\_h \frac{(\boldsymbol{u}\_i^{n+1} + \boldsymbol{u}\_i^n)}{2(\Delta\mathbf{x})^2} \quad - \frac{3}{2}\Delta\_h \frac{(\boldsymbol{F}\_i^{n+1} + \boldsymbol{F}\_i^n)}{2(\Delta\mathbf{x})^2} - \frac{1}{3}\Delta\_h^2 \frac{(\boldsymbol{u}\_i^{n+1} + \boldsymbol{u}\_i^n)}{2(\Delta\mathbf{x})^4} = \boldsymbol{f}\_i \tag{4}$$
 
$$\boldsymbol{u}\_i^0 \quad = \boldsymbol{u}(\mathbf{x}\_i, \mathbf{0})$$
 
$$\boldsymbol{D}^-(\boldsymbol{u}\_0) \quad = \quad \text{0 and } \boldsymbol{D}^+\boldsymbol{D}^-(\boldsymbol{u}\_{-1}) = \mathbf{0}$$
 
$$\boldsymbol{D}^+(\boldsymbol{u}\_N) \quad = \quad \text{0 and } \boldsymbol{D}^+\boldsymbol{D}^-(\boldsymbol{u}\_{N+1}) = \mathbf{0},$$

for *i* = 0, 1, ··· , *N*, *n* = 0, 1, ··· , *Nt*. The standard finite difference is used for forward *<sup>D</sup>*+*uni* = *<sup>u</sup>ni*+<sup>1</sup> − *uni* , backward *<sup>D</sup>*<sup>−</sup>*uni* = *uni* − *<sup>u</sup>ni*−1, central *<sup>D</sup>*0*uni* = *<sup>u</sup>ni*+<sup>1</sup> − *<sup>u</sup>ni*−1, Δ*h* = *<sup>D</sup>*+*D*<sup>−</sup>*uni* = *<sup>u</sup>ni*+<sup>1</sup> − <sup>2</sup>*uni* + *<sup>u</sup>ni*−1, and Δ2*h* = <sup>Δ</sup>*h*Δ*h<sup>u</sup>ni* = *<sup>u</sup>ni*+<sup>2</sup> − <sup>4</sup>*uni*+<sup>1</sup> + <sup>6</sup>*uni* − <sup>4</sup>*uni*−<sup>1</sup> + *<sup>u</sup>ni*−2.

Note that the nonlinear term *F* = *u*2 is linearized by taking a Taylor expansion for the nonlinear convection term *F<sup>n</sup>*+<sup>1</sup> *i*as

$$F\_i^{n+1} = F\_i^n + \Delta t \left(\frac{\partial F}{\partial t}\right)\_i^n + O(\Delta t^2),\tag{5}$$

where *∂F∂t* = *∂F<sup>∂</sup>u*Δ*<sup>u</sup>* = <sup>2</sup>*uni* (*un*+<sup>1</sup> *i* − *uni* ). Therefore, the nonlinear convection term can be linearized in the following form

$$F\_i^{n+1} + F\_i^n = 2
u\_i^n 
u\_i^{n+1}.\tag{6}$$

Substituting (6) into (4) leads to the following system

$$a\_i^u u\_{i-2}^{n+1} + b\_i^u u\_{i-1}^{n+1} + c\_i^u u\_i^{n+1} + c\_i^u u\_{i+1}^{n+1} + a\_i^u u\_{i+2}^{n+1} = d\_i^u. \tag{7}$$

where

$$a\_i^n = r, \quad b\_i^n = 2p\lambda - 3pu\_{i-1}^n - 4r - 1, \quad c\_i^n = 2p\lambda - 3p\eta\_{i+1}^n - 4r + 1 \quad c\_i^n = -4p\lambda + 6pu\_i^n + 6r$$

and

$$d\_i^{\underline{n}} = -r u\_{i-2}^{\underline{n}} + (-1 - 2p\lambda + 4r)\eta\_{i-1}^{\underline{n}} + (4p\lambda - 6r)u\_i^{\underline{n}} + (1 - 2p\lambda + 4r)u\_{i+1}^{\underline{n}} - r\_r u\_{i+2}^2 + 2\triangle x \triangle t \bar{f}\_{i+1}$$

where *p* = *t*/*<sup>x</sup>*, *r* = −*t*/3(*x*)3. This system is solved at each time step. The computational domain is taken to be b= 50 with Δ*x* = 0.1 and Δ*t* = 0.05. Our scheme is based on the linearized implicit scheme as proposed in [23,24]. This scheme is unconditionally linearly stable with a local truncation error, *<sup>O</sup>*(<sup>Δ</sup>*x*2, <sup>Δ</sup>*t*<sup>2</sup>). Therefore, there is no restriction on Δ*x* and Δ*t* (more details on this numerical scheme are found in [23,24]).
