# Analysis of the Influence of Broad Deaths on Honeybee Population
## 5 Computational Simulations
### Implicit Numerical Methods  
The choice of a numerical method is very important for the accuracy and computational cost of the solution. Usually, for complex nonlinear systems, the implicit methods perform better in terms of computational resources if a predefined accuracy is required.  
The problem with the classical implicit methods is that they are very computationally demanding, since, at each step, for each stage, we need to solve a system of nonlinear algebraic equations. It is done by a suitable iteration method, which makes the entire approach computationally intensive. A possible way to try to diminish the load is to apply a special choice of the coefficients in Equation (14). We employ the Mono-Implicit Runge-Kutta (MIRK) method [29].  
If the ODE to solve has the form  
\[\frac{\mathrm{d}y}{\mathrm{d}t}=f\big{(}t,y(t)\big{)},\quad y(0)=y_{0},\]  
then the MIRK method could be defined as follows:  
\[y_{n+1}=y_{n}+h_{n}\sum_{r=1}^{s}b_{r}K_{n+1}^{r}, \tag{14}\]  
where \(s\) is the number of stages, \(h_{n}=t_{n+1}-t_{n}\) is the size of the step, and  
\[K_{n+1}^{r}=f\big{(}t_{n}+c_{r}h_{n},(1-v_{r})y_{n}+v_{r}y_{n+1}+h_{n}\sum_{j= 1}^{r-1}a_{rj}K_{n+1}^{j}\big{)}.\]  
The other coefficients could be defined by means of the special form of the Butcher tableau in Equation (15):  
\[\begin{array}{c|ccccccccc}c_{1}&v_{1}&0&0&\cdots&0&0\\ c_{2}&v_{2}&a_{2,1}&0&\cdots&0&0\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ c_{s}&v_{s}&a_{s,1}&a_{s,2}&\cdots&a_{s,s-1}&0\\ \hline&&b_{1}&b_{2}&\cdots&\cdots&b_{s}\end{array}\equiv\begin{array}{c|c |c}C&V&A\\ \hline&B\\ \hline&B\\ \end{array}. \tag{15}\]  
For consistency, it must be true that \(C=V+A*e\), where the elements of vector \(C=(c_{1},c_{2},\ldots,c_{s})^{\top}\) are called abscissae, \(B=(b_{1},b_{2},\ldots,b_{s})\), \(\sum_{i=1}^{s}b_{i}=1\), \(V=(v_{1},v_{2},\ldots,v_{s})^{\top}\), \(e=(1,1,\ldots,1)^{\top}\), and \(A=\{a_{ij}\}_{i,j=1}^{s}\) is the Runge-Kutta matrix. The abscissae fulfill \(c_{r}=v_{r}+\sum_{i=1}^{s}a_{ri}\).  
The benefit of the MIRK over the other implicit methods is that the stages explicitly relate to each other, and they appear implicit only with respect to \(y_{n+1}\). In this manner, the method retains the high accuracy of the implicit class, while the computational resources needed are significantly reduced.  
In this case, our DDE looks like  
\[\frac{\mathrm{d}y}{\mathrm{d}t}=f\big{(}t,y(t),y(t-\tau)\big{)},\quad y(t)=y_ {0}(t)\text{ for }-\tau\leq t\leq 0,\]  
then it is possible to adapt the MIRK method for delayed equations [30].  
Similar to Equation (14), the delayed MIRK is defined by  
\[y_{n+1}=y_{n}+h_{n}\sum_{r=1}^{s}b_{r}K_{n+1}^{r},\]  
but now the stages are defined as  
\[K_{n+1}^{r}=f\big{(}t_{n}+c_{r}h_{n},(1-v_{r})y_{n}+v_{r}y_{n+1}+h_{n}\sum_{j=1} ^{r-1}a_{rj}K_{n+1}^{j},y(t_{n}+c_{r}h_{n}-\tau)\big{)}.\]  
This modification is enough if the delay \(\tau\) is multiple of the step \(h_{n}\), since, then, the delayed term \(y(t_{n}+c_{r}h_{n}-\tau)\) would always coincide with a grid point.  
It is possible to obtain different schemes, if the coefficients \(A\), \(V\), and \(B\) vary. Of course, they need to meet certain conditions [31]. In our experiments, we used the third-order scheme [30].  
\[\begin{array}{c|c|c|c}0&0&0&0&0&0\\ 1&1&0&0&0&0\\ \frac{3}{4}&\frac{27}{32}&\frac{3}{64}&-\frac{9}{64}&0&0\\ \frac{1}{5}&\frac{13}{125}&\frac{16}{125}&-\frac{4}{125}&0&0\\ \hline&&\frac{5}{9}&-\frac{3}{8}&\frac{128}{99}&-\frac{125}{264}\end{array}.\]