# Numerical Simulation of the Kelvin Wake Patterns
## 3 Numerical Methods
### Jacobian\(-\)Free Newton\(-\)Krylov Method  
After numerical discretization, a nonlinear system of equations could be obtained, as follows:  
\[\mathbf{F}(\mathbf{u})=0, \tag{35}\]  
where \(\mathbf{u}\) is the vector of unknowns of the length \(2(N+1)M\). This nonlinear system can be solved by [PERSON]'s method [35]. The iterative process is as follows: set an initial guess \(\mathbf{u_{0}}\), for \(t=0,\ 1,\ 2,\ \ldots\), until convergence:  
\[\mathbf{J}(\mathbf{u_{t}})\mathbf{\delta u_{t}}=-\mathbf{F}(\mathbf{u_{t}}), \tag{36}\]  
\[\mathbf{u_{t+1}}=\mathbf{u_{t}}+\mathbf{\delta u_{t}}, \tag{37}\]  
where \(\mathbf{J}(\mathbf{u_{t}})=\mathbf{F}(\mathbf{u_{t}})\) is the system Jacobian [38]. However, [PERSON]'s method cannot solve large nonlinear systems efficiently. Therefore, based on [PERSON]'s method, this paper adopts the Jacobian\(-\)free Newton-Krylov method with a preconditioner, which is carefully constructed with entries taken from the Jacobian of the linearized problem. There is one more inner iteration, since Equation (36) is solved by using the iterative Generalized Minimum Residual algorithm [39]. The GMRES method is one of the Krylov subspace methods that are attractive as linear solvers in the context of nonlinear Newton iteration. Since each iteration does not require the exact value of \(\mathbf{\delta u_{t}}\) in solving nonlinear equations, the GMRES method can increase computation speed. Firstly, the approximate solution of \(\mathbf{\delta u_{t}}\) is found by projecting obliquely onto the Krylov subspace  
\[\mathcal{K}_{m}\Big{(}\mathbf{J_{t}P^{-1},F_{t}}\Big{)}=\ \text{span}\{ \mathbf{F_{t}},\mathbf{J_{t}P^{-1}F_{t}},\ \ldots,\ \Big{(}\mathbf{J_{t}P^{-1}}\Big{)}^{m-1}\mathbf{F_{t}}\} \tag{38}\]  
where \(m\) is the value of the subspace dimension and the accuracy of the solution increases with the subspace dimension, \(\mathbf{J_{t}}=\mathbf{J}(\mathbf{u_{t}}),\mathbf{F_{t}}=\mathbf{F}(\mathbf{u_{t}})\). The matrix \(\mathbf{P}\approx\mathbf{J_{t}}\) is the preconditioner matrix, whose purpose is to construct an approximation to the Jacobian \(\mathbf{J_{t}}\), which is cheap to form and to factorize. The calculation speed of the GMRES method can be significantly improved with a preconditioner matrix, because the spectrum of the preconditioned Jacobian \(\mathbf{J_{t}P^{-1}}\) exhibits a clustering of eigenvalues [40]. An initial linear residual \(\mathbf{r_{0}}\) is defined, given an initial guess \(\mathbf{\delta u_{0}}\), for the Newton correction,  
\[\mathbf{r_{0}}=-\mathbf{F}(\mathbf{u_{0}})-\mathbf{J_{0}P^{-1}\delta u_{0}}. \tag{39}\]Subsequently, the GMRES iteration minimizes \(\left\|\mathbf{r_{t}}\right\|_{2}\) to a suitable value. These Jacobian\(--\) vector products can be approximated by applying first\(--\)order difference quotients:  
\[\mathbf{J}_{t}\mathbf{P}^{-1}\mathbf{v}\approx\frac{\mathbf{F}\big{(}\mathbf{u_{t}}+h \mathbf{P}^{-1}\mathbf{v}\big{)}-\mathbf{F}(\mathbf{u_{t}})}{h} \tag{40}\]  
where \(\mathbf{v}\) represents an arbitrary vector used in building the Krylov subspace [41], and the \(h\) is a small perturbation [40]  
\[h=\frac{\sqrt{(1+\left\|\mathbf{u_{t}}\right\|)h_{mach}}}{\left\|\mathbf{v}\right\|}. \tag{41}\]  
In order to calculate Equation (36) conveniently, the preconditioner can be divided up into four equal submatrices and factorized using the block decomposition,  
\[\mathbf{P}= \begin{bmatrix}\mathbf{A}&\mathbf{B}\\ \mathbf{C}&\mathbf{D}\end{bmatrix} = \begin{bmatrix}\mathbf{I}&0\\ \mathbf{C}\mathbf{A}^{-1}&\mathbf{I}\end{bmatrix}\begin{bmatrix}\mathbf{A}&0 \\ 0&\mathbf{D}-\mathbf{C}\mathbf{A}^{-1}\mathbf{B}\end{bmatrix}\begin{bmatrix} \mathbf{I}&\mathbf{A}^{-1}\mathbf{B}\\ 0&\mathbf{I}\end{bmatrix} , \tag{42}\]  
where \(\mathbf{I}\) is the unit matrix, \(\mathbf{A},\ \mathbf{B},\ \mathbf{C}\) and \(\mathbf{D}\) are the four submatrices, which are constructed in the base of Jacobian \(\mathbf{J}_{t}\). It is possible to solve the system \(\mathbf{P}\mathbf{r}=\mathbf{v}\) by performing the following operations,  
\[