Referring to Figure 64, we use the following numbering: $$ \begin{align} \tag{4.153} \begin{array}{l} x_i = (i-1)\;h \ , \ i=1,2,...,N+1 \textrm{ with } h = \frac{1}{N}\\ M = \frac{N}{2}+1 \textrm{ with N being an even number and }x_M = 0.5 \end{array} \end{align}$$

With \( k = k_a \) for \( i = M \) we get that \( \beta = \beta_a \) for \( i=1,2,...,M \) and \( \beta = \beta_b \) for \( i > M \). For \( i\not= M \) we can write (4.150) $$ \begin{equation} \tag{4.154} \frac{d}{dx}\left( x\frac{d\theta}{dx}\right) - \beta^2 \;\theta(x) = 0 \end{equation} $$

Discretizing (4.154) using (2.45) in Chapter (2 Initial value problems for Ordinary Differential Equations): $$ \begin{equation*} -x_{i- \frac{1}{2}} \;\theta_{i-1} + (x_{i+ \frac{1}{2}} + x_{i- \frac{1}{2}} + \beta^2 h^2) \;\theta_i - x_{i + \frac{1}{2}} \;\theta_{i+1} = 0 \end{equation*} $$

Noting that \( x_{i- \frac{1}{2}} = (i - \frac{3}{2}) \;h \) and \( x_{i+ \frac{1}{2}} = (i - \frac{1}{2})\;h \) we can write: $$ \begin{equation*} \left(i - \frac{3}{2}\right) \;\theta_{i-1} + [2(i-1) + \beta^2h] \;\theta_i - \left(i- \frac{1}{2} \right) \;\theta_{i+1} = 0 \end{equation*} $$

or divided by \( i \): $$ \begin{equation} \tag{4.155} \ - \left( 1-\frac{3}{2i} \right) \;\theta_{i-1} + \left[ 2 \;\left( 1 - \frac{1}{i} \right) + \frac{\beta^2h}{i} \right] \;\theta_i - \left(1 - \frac{1}{2i} \right) \;\theta_{i+1} = 0 \end{equation} $$

(4.155) is used for \( i = 2,3,...,N \), however, except for \( i = M = N/2+1 \).

For i \( = \) 1

From (4.82) we get the boundary condition for \( x=0 \): $$ \begin{equation} \dfrac{d\theta}{dx}(0)= \beta_a^2 \;\theta(0) \to \dfrac{-3 \theta_1 + 4 \theta_2 - \theta_3}{2h}= \beta_a^2 \;\theta_1, \tag{4.156} \end{equation} $$

where we have used forward differences (See forward).

This results in: $$ \begin{equation} \tag{4.157} \theta_3 = 4\theta_2 - (3 + 2h\beta_a^2) \;\theta_1 \end{equation} $$

Writing (4.155) for \( i= 2 \): $$ \begin{equation*} -\left(1- \frac{3}{4}\right) \;\theta_1 + \left[ 2\;\left(1 - \frac{1}{2}\right) + \frac{\beta_a^2h}{2}\right] \;\theta_2 - \left(1 - \frac{1}{4}\right) \;\theta_3 \end{equation*} $$

and considering (4.151) and (4.157) results in: $$ \begin{equation} \tag{4.158} \left( 1 + \frac{3h}{2}\right) \;\theta_1 - \left( 1 - \frac{h}{2}\right) \;\theta_2 = 0 \end{equation} $$

This is the first equation of our system.

For i \( = \) M.

According to the discretization illustrated in Figure 65, here we use an equation of the form given in (4.150), which discretized using (2.45) from Chapter (2 Initial value problems for Ordinary Differential Equations), results in: $$ \begin{align*} \begin{array}{c} \beta_{M- \frac{1}{2}}^2 = \beta_a^2 = 2.0 \ , \ \beta_{M+\frac{1}{2}}^2 = \beta_b^2 = 8.0\\ x_{M-\frac{1}{2}} = \dfrac{h}{2}(N-1)= \dfrac{(1-h)}{2} \ , \ x_{M+\frac{1}{2}} = \dfrac{h}{2}(N+1) = \dfrac{(1+h)}{2} \end{array} \end{align*} $$

which provides: $$ \begin{equation} \tag{4.159} -4 \;(1-h) \;\theta _{M-1} + [5-3h+16h^2]\;\theta_M - (1+h) \;\theta_{M+1} = 0 \end{equation} $$

For i \( = \) N we can use (4.155) with \( \beta^2 = \beta_b^2 = 8.0 \): $$ \begin{equation} \tag{4.160} -\left(1-\frac{3h}{2}\right) \;\theta_{N-1} + 2 \;[(1-h)+4h^2] \;\theta_N = \left(1-\frac{h}{2}\right) \end{equation} $$

Equations (4.158)–(4.160) constitute a linear system of equations with a tridiagonal matrix and can be solved as such using the Thomas algorithm. In addition to the temperature, we also want to calculate the temperature gradient \( \theta'(x) \) . For \( x=0 \) we use (4.82) and interface condition \( \frac{d\theta}{dx}(0.5+)\equiv \left( \frac{d\theta}{dx}\right)_b \) found from (4.149). The other values are computed using standard central differences, besides for x \( = \) 0.5 og x \( = \) 1.0, where second order backward differences are used. Let us look at the use of the differential equation as an alternative.

We integrate (4.154): $$ \begin{equation} \tag{4.161} \frac{d\theta_i}{dx}=\frac{\beta^2}{x_i} \int_{x_1}^{x_i} \theta(x)dx \ , \ i= 2,3,...,N+1 \end{equation} $$

Integral \( I = \int_{x_1}^{x_i}\theta(x)dx \) can be calculated using the trapezium method, for example . The computation of (4.161) is shown below in pseudo-code form: $$ \begin{align*} \begin{array}{rrrlll} &&\theta_1' &:=& \beta_a^2\;\theta_1\\ &&s&:=& 0\\ &&\text{Utfør for } i&:=&2,...,N+1\\ &&x&:=&h \; (i-1)\\ &&s&:=&s+0.5h \; (\theta_i+\theta_{i-1})\\ \end{array} \end{align*} $$ $$ \begin{equation*} \text{Dersom $(i\le M)$ sett }\theta_i ':= \frac{\beta_a^2\;s}{x} \text{ ellers }\theta_i':= \frac{\beta_b^2 \;s}{x} \end{equation*} $$

In the following table we have used (4.161) and the trapezium method. In this case, the accuracy of the two methods for calculating \( \theta'(x) \) is quite similar as both \( \theta \) and \( \theta' \) are smooth functions (except at x \( = \) 0.5), but generally the integration method will be more accurate if there are discontinuities..

Solution of equation (4.150) with \( h = 0.01 \)