## 2 The Examined Numerical Methods
### Implementations of the Convection Term  
Until this point, the LH algorithm has been given only for the conduction and the heat source term. Now, the \(-Ku\) term is also included, which is done in several ways.  
1. The explicit treatment means that one calculates the increment due to the term and simply adds it to the final value of the new \(u\), just as it would happen for the explicit Euler method. For example, in the case of the first stage with a halved time step, the increment is \(-\frac{\Delta t}{2}Ku_{i}^{0}\), thus we have  
\[u_{i}^{1/2}=\frac{u_{i}^{0}+\
icefrac{{r}}{{2}}\big{(}u_{i-1}^{0}+u_{i+1}^{0} \big{)}+\
icefrac{{\Delta t}}{{2}}\cdot\
icefrac{{q_{i}}}{{2}}}{1+r}- \
icefrac{{h}}{{2}}Ku_{i}^{0}. \tag{8}\]  
We note that we exemplify the versions with a first-stage formula, since it is the least nontrivial due to the half-sized time step.  
2. Quasi-exact treatment means that we analytically solve the ODE  
\[\frac{du}{dt}=-K\cdot u, \tag{9}\]  
and then take the effect of the convection terms into account in a separate calculation at the end of each stage. For example, in the case of the first stage with a halved time step, we have  
\[u_{i}^{1/2,\text{ temp}}=\frac{u_{i}^{0}+\
icefrac{{r}}{{2}}\big{(}u_{i-1}^{0}+u_{i +1}^{0}\big{)}+\
icefrac{{h}}{{2}}\cdot\
icefrac{{q_{i}}}{{2}}}{1+r} \tag{10}\]  
and  
\[u_{i}^{1/2}=e^{-Kh/2}u_{i}^{1/2,\text{ temp}}. \tag{11}\]  
We expect that this version has an outstanding stability, since the absolute value of the solution is always smaller than the temporary value \(u_{i}^{1/2,\text{ temp}}\), where only the conduction and the source terms are taken into account.  
3. Pseudo-implicit treatment means that the \(u\) variable in the convection term is taken into account at the new time level, so Equation (4) is extended as follows:  
\[u_{i}^{n+1}=u_{i}^{n}+r\Big{[}u_{i-1}^{n}+u_{i+1}^{n}-2\theta u_{i}^{n}-2(1- \theta)u_{i}^{n+1}\Big{]}-Khu_{i}^{n+1}. \tag{12}\]  
With this, the \(K\) term turns up only in the denominator, thus the first stage formula is the following:  
\[u_{i}^{1/2}=\frac{u_{i}^{0}+\
icefrac{{r}}{{2}}\big{(}u_{i-1}^{0}+u_{i+1}^{0} \big{)}+\
icefrac{{h}}{{2}}\cdot\
icefrac{{q_{i}}}{{2}}}{1+r+hK/2}. \tag{13}\]  
Note that in our last paper [31], only this implementation was proposed.  
4. Now, the \(u\) variable in the convection term is taken into account at the old time level, so in Equation (12), the last term is changed to \(-Khu_{i}^{n}\) This means that the \(K\) term turns up only _inside_ the numerator, so the first stage formula is the following:  
\[u_{i}^{\sfrac{1}{2}}=\frac{u_{i}^{0}+\
icefrac{{r}}{{2}}\big{(}u_{i-1}^{0}+u_{i+ 1}^{0}\big{)}+\
icefrac{{h}}{{2}}\cdot q_{i}-Khu_{i}^{0}/2}{1+r}. \tag{14}\]  
Due to the lack of a better name, we call this version temporarily 'inside'.  
5. Mixed treatment means that we make a linear combination of the last two versions (pseudo-implicit and inside) at the level of Equation (12), where the last term is changed to \(-pKhu_{i}^{n+1}-(1-p)Khu_{i}^{n}\). The real parameter \(p\) has a similar role as that which \(\theta\) has for the conduction case. Now, the formulas are the following:  
\[\begin{split}\text{First stage}:&\;u_{i}^{\sfrac{1 }{2}}=\frac{u_{i}^{0}+\
icefrac{{r}}{{2}}\big{(}u_{i-1}^{0}+u_{i+1}^{0}\big{)}+ \
icefrac{{h}}{{2}}\cdot q_{i}-(1-p)Khu_{i}^{0}/2}{1+r+pKh/2},\\ \text{Intermediate stages}:&\;u_{i}^{n+1}=\frac{(1-r)u_{ i}^{n}+\big{(}u_{i-1}^{n+1}+u_{i+1}^{n+1}\big{)}+\
icefrac{{h}}{{2}}\cdot q_{i}-(1-p)Khu_{i}^{n}}{1+r+ pKh},\\ \text{Last stage}:&\;u_{i}^{T}=\frac{(1-\
icefrac{{r }}{{2}})u_{i}^{T-\
icefrac{{r}}{{2}}}+\
icefrac{{r}}{{2}}\big{(}u_{i-1}^{T}+u_ {i+1}^{T}\big{)}+\
icefrac{{h}}{{2}}\cdot q_{i}-(1-p)Khu_{i}^{T-\
icefrac{{r }}{{2}}}/2}{1+\
icefrac{{r}}{{2}}+pKh/2}.\end{split} \tag{15}\]  
We performed tests with several values of \(p\), but here we present them only for three values, namely, \(p=\
icefrac{{1}}{{3}}\), \(\
icefrac{{1}}{{2}}\), \(\
icefrac{{2}}{{3}}\).