## 4 Results  
### 1D Case  
#### 4.1.1 Convergence Analysis  
The theoretical order of accuracy of the numerical schemes that we use here is _first_ -order in time and _second_ order in space for the backward Euler and _second_ order in time and _second_ order in space for the Crank-Nicolson method. Of course, these theoretical orders exist for the linear models and linear methods. Whenever the model becomes nonlinear or sharp structures appear in the solution, it may become difficult to verify the expected order of accuracy. We have both situations in our model and solution profiles. Therefore, it makes sense to check whether we lose numerical accuracy in time or space. To do this, we will perform temporal and spatial convergence analysis, that subsequently described. Although there exists analytical solutions to the Fisher-Kolmogorov equation in certain situations [44], we prefer computational means to carry the convergence analysis. This is mainly because our numerical analysis can be easily applied to more general cases such as a coupled nonlinear fluid dynamical system with more physical effects. The time convergence analysis was performed in the following manner. We fix the number of grid points (e.g., \(M=2048\)) and run the code with sequentially refined time steps up to the final time \(t_{final}=10\). For instance, we obtain the following approximate solution set \(c^{\Delta t},c^{\Delta t},c^{\Delta t},\ldots,c^{\Delta t}\) where \(k=0,\,1,\ldots,\,8\) and \(\Delta t=10^{-1}\). Then, we measure the \(L_{p}\) norms of errors between two consecutive time step solutions, e.g.,  
\[||c^{\Delta t}\,-\,c^{\Delta t}||_{p},\,||c^{\Delta t}\,-\,c^{\Delta t}||_{p},\ldots,\]  
where \(p\in\{1,\,2,\,\infty\}\). These time errors can be superimposed (in log scale) on top of the reference first- and second-order slope indicator lines to provide insights into the convergence behavior. Alternatively, one can define the time rate of convergence by  
\[r=\ln\Bigl{(}||E(\Delta t/2^{k-1})||/||E(\Delta t/2^{k})||\Bigr{)}\,/\ln 2, \tag{25}\]  
where  
\[E(\Delta t/2^{k-1})=||c^{\Delta t/2^{k-1}}\,-\,c^{\Delta t/2^{k}}||_{p}\,\, \,\,\,\,\,\,\,k=1,\ldots,\,8.\]  
Based on the second description, we made convergence tables to demonstrate and verify the theoretical time order of accuracy of the numerical schemes. Table 1 illustrates the time behavior of the Crank-Nicolson method using different norms. The convergence rates slightly change depending on the norm, however, they are all very close to the theoretical value which is _two_. In Table 2, we verify the time convergence behavior of the backward Euler method. Naturally, time errors are bigger compared to the Crank-Nicolson results, but again, we have a perfect match between the theoretical and numerical rate of convergence. This verification analysis is very important because it not only means that we have a reliable code implementation, but it also numerically proves that the JFNK method consistently converges all nonlinearities in the physical system as intended. To numerically demonstrate the spatial convergence, we perform _nine_ runs using \(\frac{\Delta x}{2^{k}}\)\((\,k\,=\,0,\,1,\ldots,\,8)\) and a fixed fine time step \(\Delta t\,=\,10^{-5}\) where \(\Delta x=\frac{50}{2^{6}}\). In these runs, the time scheme is set to be the CN method that, by construction, contributes smaller time errors to the truncation term. With these settings and similar error norm measures described above, we formed Table 3 at the final time \(t_{final}\,=\,10\). Table 3 clearly shows second-order spatial convergence. Again, this is a good indication of reliable and accurate code implementation. This also means that we can readily extend our work to multi-dimensional and more complex models. Here, we note that one could fix \(\Delta t/\Delta x^{2}\) or \(\Delta t^{2}/\Delta x^{2}\) and perform only one convergence analysis (space-time) for each numerical scheme to check the possible accuracy loss.