# A study of the One-Dimensional
Heat-Conduction Equation with Radiation

Mihai Halic

###### Abstract.

We consider a boundary value problem (BVP) modelling one-dimensional heat-
conduction with radiation, which is derived from the Stefan-Boltzmann law. The
problem strongly depends on the parameters, making difficult to estimate the
solution. In here we apply an analytical approach to determine upper and lower
bounds to the exact solution of the BVP, which allows estimating the latter.
Finally, we support our theoretical arguments with numerical data, by
implementing them into the MAPLE computer program.

###### Key words and phrases:

heat-conduction, radiation, boundary layer, envelopes, Runge-Kutta

###### 2010 Mathematics Subject Classification:

Primary 34B15; Secondary 34L30, 34B60

##  1\. Introduction

Firnett-Troesch [4] applied the so-called ‘shooting method’ to numerically
analyse two second order BVPs, which depend in a sensitive way on the
parameters. The author devised [5] an analytic approach for determining
explicit approximate solutions to such BVPs, and compared the outcome with
existing numerical results for the ‘first Troesch-equation’
$y^{\prime\prime}=L\sinh(Ly)$. The analytical method turned out to be several
orders of magnitude preciser than previous numerical methods, especially for
large values of $L$. Here we develop this approach, and consider the ‘second
Troesch-equation’ that is,

$\hskip 43.36464ptu^{\prime\prime}_{\rm heat}=b^{2}(u^{4}_{\rm
heat}-t^{4}),\quad u_{\rm heat}(0)=1,\;u_{\rm heat}(1)=t\in(0,1),\hskip
43.36464pt\text{(HCR)}$

which describes the steady-state, one-dimensional heat conduction with
radiation. The equation depends on two parameters: $b$ is related to the
radiative properties of materials, while $t$ is roughly the ratio between the
initial and final temperatures.

Our goal here is to approximately solve (HCR), understand the dependence of
its solution on the parameters, elucidate the law governing the boundary
layer, etc. These objectives turn out to be challenging. There seems to be a
rather restricted literature dedicated to the analysis of this very equation,
in spite of being interesting:

  * –

On the mathematics side, the ODE is far from being classically integrable.
Existence results are proved in [8, 7, 10]. Power series expansions are
obtained in [11, 12]; their terms are complicated and, for good
approximations, the order of the development quickly increases with the
parameters.

From a numerical perspective, (HCR) often has the so-called ‘boundary layer
property’, $u_{\rm heat}$ has huge derivative near $x=0$, which is a major
obstacle. For running the shooting method, one needs a very good estimate of
either $u^{\prime}_{\rm heat}(0)\ll 0$ or $u^{\prime}_{\rm heat}(1)\approx 0$.
Unfortunately, for larger $b$ —physically realistic situations involve
$b\approx 30\dots 40$— these values are hard to estimate. Indeed, Firnett-
Troesch consider only a few values ($\,0.8,0.4,0\,$) of $t$, and deal with
values of $b$ up to $20$.

[4, p.427]: |  [4, Figure 6.1, p.426]:
---|---
“If the differential equation is integrated with a guess […] which is even slightly in error, then a strongly diverging solution is obtained.” |

Thus one faces a chaotic situation —tiny changes lead to major variations—,
which raises the issue how to deal with a non-integrable ODE, without
completely relying on computer power. Our approach is to determine explicit
upper and lower bounds —so-called envelopes— and deduce the dependence on the
parameters. Quite remarkably, the conclusions follow naturally.

  * –

On the physics side, the equation is based on the Stefan-Boltzmann law for the
heat radiation; it’s a rare physical law involving $4^{th}$-power dependence
on a parameter (the temperature). Numerous textbooks [1, 2, 3, 6, 9], written
from both theoretical and engineering standpoints, show the importance of heat
conduction and radiation.

### Outline

At the core of our approach, there is a careful analysis of the heat
conduction equation. The large variations of the output, generated by
insignificant changes in the input, can’t be controlled without taming the
mechanism which yields solution. The computer work supplements and enhances
the presentation, also reinforces the (possibly confusing) theoretical
arguments with numerical evidence.

The article is organized as follows: we start by estimating the solution of a
transcendental equation, which occurs in our analysis. In Section 3 we obtain
the first estimates for $u_{\rm heat}$, actually an upper bound. This estimate
is already crucial:

  * –

It yields an approximate value for the initial derivative, valid for any
choices of the parameters $b$ and $t$, which is otherwise impossible to guess.
This approximate value allows running the Runge-Kutta method for (HCR), with
quite a remarkable precision, even for low values of the parameters.111The
author reported on the application of the Runge-Kutta method, briefly
discussed here in Section 3.4, at the ICNAAM 2023 conference. The short
abstract of the report will appear in the AIP Conference Proceedings.

  * –

It allows pinning down the law governing the boundary layer property, in terms
of the parameters $b,t$. Remarkably, this law has a geometric and a
variational interpretation.

In the Sections 4, 5, 6, we refine our considerations, and determine accurate
envelopes for the exact solution of the heat equation. Finally, we implement
our considerations into the MAPLE software, in Section 7; the code is given,
the interested reader may freely experiment with the favourite values of the
parameters. Numerical data shows that the precision of the estimates increases
with the slope of the boundary layer.

  * Acknowledgments

The author thanks K. Mustapha for discussions related to numerical modelling,
which helped improving the rule governing the boundary layer.

##  2\. Background matter

The following calculations are tedious, but essential.

###  2.1. A transcendental equation I

We consider real numbers $\zeta>1$, ${\Lambda}>0$, and the equation in the
unknown $c\geqslant 0$:

  * $\displaystyle\begin{array}[]{r}\displaystyle\frac{\zeta-c}{\zeta+c}=\frac{1-c}{1+c}\cdot e^{2{\Lambda}c}\;\Leftrightarrow\;(\zeta-1)\cdot\frac{c}{\zeta-c^{2}}=\tanh({\Lambda}c)\;\Leftrightarrow\;\zeta=\underbrace{c\cdot\frac{1-c\cdot\tanh({\Lambda}c)}{c-\tanh({\Lambda}c)}}_{rhs(c)}.\end{array}$ (2.1)

First we analyse the poles of the right-hand side. The equation
$c=\tanh({\Lambda}c)$ has the solution $c_{0}({\Lambda})=0$, for
${\Lambda}\in(0,1]$. For ${\Lambda}>1$, it admits a (unique) solution
$c_{0}({\Lambda})\in(0,1)$. In all cases, $c_{0}({\Lambda})$ can be expressed
as an infinite iteration of the ‘tanh’ function:

$c_{0}({\Lambda})=\tanh({\Lambda}\cdot\tanh({\Lambda}\cdot\tanh({\Lambda}\cdot\dots))).$

This formula allows computing approximate values for $c_{0}({\Lambda})$:

  * $\displaystyle\begin{array}[]{|l|l|}\hline\cr{\Lambda}&c_{0}({\Lambda})=\tanh({\Lambda}\cdot\tanh({\Lambda}\cdot\tanh({\Lambda}\cdot\dots)))\\\ \hline\cr{\Lambda}=1&c_{0}(1)=0\\\ \hline\cr{\Lambda}=1.10&c_{0}(1.10)\approx 0.5\\\ \hline\cr{\Lambda}=1.50&c_{0}(1.50)\approx 0.85\\\ \hline\cr{\Lambda}=1.67&c_{0}(1.67)\approx 0.9\\\ \hline\cr{\Lambda}=2&c_{0}(2)\approx 0.95\\\ \hline\cr{\Lambda}\geqslant 5&c_{0}({\Lambda})>0.999\approx 1\\\ \hline\cr\end{array}$ (2.2)

A numerical experiment shows that $20$ iterations of the function
$\tanh({\Lambda}\cdot)$ are necessary to obtain a two digit approximation —an
upper bound— of $c_{0}({\Lambda})$, for ${\Lambda}\geqslant 1.1$.

  * Lemma 2.1

Given $\zeta>1$ and ${\Lambda}>0$, the equation (2.1) admits a unique solution
$c=c(\zeta,{\Lambda})$, as follows:

    1. (i)

for ${\Lambda}>1$, one has $c\in(c_{0}({\Lambda}),1)$.

    2. (ii)

for ${\Lambda}\in(0,1]$, $\zeta<\frac{1}{1-{\Lambda}}$, one has $c\in(0,1)$.
(For ${\Lambda}=1$, we let $1/(1-1)=+\infty$.)

###### Proof.

In the first case, we have $\mathop{\lim}\limits_{c\to
c_{0}({\Lambda})^{+}}rhs(c)=+\infty$ and $rhs(1)=1<\zeta$. In the second case,
the denominator doesn’t vanish, and we have $\mathop{\lim}\limits_{c\to
0^{+}}rhs(c)=\frac{1}{1-{\Lambda}}>\zeta$, $rhs(1)=1<\zeta$. Since $c\mapsto
rhs(c)$ is decreasing, the uniqueness follows. ∎

###  2.2. Bounding the solution c

We determine explicit upper and lower bounds for the exact solution of (2.1),
depending on the parameters $\zeta,{\Lambda}$. The previous Lemma shows that
we should distinguish between two cases: ${\Lambda}>1$ (actually, we take
${\Lambda}\geqslant 1.1$) and ${\Lambda}<1.1$. The approaches to estimating
the solutions are different.

####  2.2.1. Case ${\Lambda}\geqslant 1.5$

Polynomial approximations fail in this situation because the $\tanh$-function
approaches exponentially fast the value $1$. Since $c$ belongs to $(0,1)$, we
write:

$\;c=\tanh(Y)=\tanh({\Lambda}X)>c_{0}({\Lambda})=\tanh({\Lambda}\cdot
c_{0}({\Lambda})),\;X>c_{0}({\Lambda}).$

The addition formula for the $\tanh$-function yields:

  * $\displaystyle\zeta=\frac{c}{\tanh(Y-{\Lambda}\tanh(Y))}=\frac{\tanh({\Lambda}X)}{\tanh({\Lambda}(X-\tanh({\Lambda}X)))},\quad X=c_{0}({\Lambda})+\gamma,\;\gamma>0.$ (2.3)

To determine explicit approximations, we consider the inequalities
‘$\leqslant$’ and ‘$\geqslant$’.

  * ‘$\leqslant$’

Since $c>c_{0}({\Lambda})$, it’s enough to have
$\tanh\big{(}{\Lambda}(X-\tanh({\Lambda}X))\big{)}\leqslant\zeta^{-1}c_{0}({\Lambda})$,
so

$X-\tanh({\Lambda}X)\leqslant{\Lambda}^{-1}\cdot\tanh^{-1}(\zeta^{-1}\cdot
c_{0}({\Lambda})).$

By using (2.2.1), the left-hand side is less than
$\gamma+c_{0}({\Lambda})-\tanh({\Lambda}c_{0}({\Lambda}))=\gamma$. Overall, we
deduce that the following value is convenient:

    * $\displaystyle c_{-}(\zeta,{\Lambda})=\tanh\biggl{(}{\Lambda}\cdot c_{0}({\Lambda})+\tanh^{-1}(\zeta^{-1}\cdot c_{0}({\Lambda}))\biggr{)}.$ (2.4)

  * ‘$\geqslant$’

As $c<1$, it’s enough to have:

$\begin{array}[]{l}\tanh\big{(}{\Lambda}(X-\tanh({\Lambda}X))\big{)}\geqslant\zeta^{-1}\geqslant\zeta^{-1}c,\\\\[8.61108pt]
\gamma+\tanh({\Lambda}c_{0}({\Lambda}))-\tanh({\Lambda}X)=X-\tanh({\Lambda}X)\geqslant{\Lambda}^{-1}\tanh^{-1}(\zeta^{-1}).\end{array}$

A concavity argument yields
$\tanh({\Lambda}X)-\tanh({\Lambda}c_{0}({\Lambda}))\leqslant\frac{{\Lambda}\gamma}{\cosh^{2}({\Lambda}c_{0}({\Lambda}))}$,
which implies that the following value is convenient:

    * $\displaystyle c_{+}(\zeta,{\Lambda})=\tanh\biggl{(}{\Lambda}\cdot c_{0}({\Lambda})+\frac{\tanh^{-1}(\zeta^{-1})}{1-\frac{{\Lambda}}{\cosh^{2}({\Lambda}\cdot c_{0}({\Lambda}))}}\biggr{)}.$ (2.5)

Note that, by differentiating
$c_{0}({\Lambda})=\tanh({\Lambda}c_{0}({\Lambda}))$ for ${\Lambda}$, one finds
that the denominator above is positive.

####  2.2.2. Case ${\Lambda}\geqslant 5,\;\zeta\leqslant e^{{\Lambda}}$

In this situation, the previous formulae can be simplified, by eliminating
$c_{0}({\Lambda})\approx 1$: we replace $X=1+\gamma$ in (2.2.1), with
$\gamma\geqslant 0$. The very same approach yields:

  * $\displaystyle\begin{array}[]{rl}\text{`$\leqslant$' is satisfied by}&\displaystyle c_{-}(\zeta,{\Lambda}):=\tanh\Big{(}{\Lambda}\cdot\tanh({\Lambda})+\tanh^{-1}\big{(}\zeta^{-1}\cdot\tanh({\Lambda})\big{)}\Big{)};\\\\[8.61108pt] \text{`$\geqslant$' is satisfied by}&\displaystyle c_{+}(\zeta,{\Lambda}):=\tanh\Biggl{(}\frac{\tanh^{-1}(\zeta^{-1})+\frac{{\Lambda}(\sinh(2{\Lambda})-2{\Lambda})}{2\cosh^{2}({\Lambda})}}{1-\frac{{\Lambda}}{\cosh^{2}({\Lambda})}}\Biggr{)}.\end{array}$ (2.6)

The condition $\zeta\leqslant e^{\Lambda}$ is necessary to ensure that
$\gamma$ is positive.

####  2.2.3. Case ${\Lambda}\in(0,1.5)$

In this case, the bounds $c_{0}({\Lambda})<c<1$ are imprecise, we need better
ones. The series development of $\frac{\tanh(x)}{x}$ yields the estimate:

  * $\displaystyle\begin{array}[]{c}\Bigl{|}\frac{\tanh(x)}{x}-(1-0.2x^{2})\Bigr{|}\leqslant 0.15x^{2},\quad\forall\,x\in[0,1.5].\end{array}$ (2.7)

Thus the solutions of the bi-quadratic equations

  * $\displaystyle\begin{array}[]{l}\frac{\zeta-1}{{\Lambda}}=(\zeta-c^{2})\cdot\Bigl{(}1-0.35\cdot{\Lambda}^{2}c^{2}\Bigr{)},\quad\frac{\zeta-1}{{\Lambda}}=(\zeta-c^{2})\cdot\Bigl{(}1-0.05\cdot{\Lambda}^{2}c^{2}\Bigr{)},\\\\[8.61108pt] \displaystyle c_{\rho}(\zeta,{\Lambda})={\Biggl{[}\frac{2\frac{1-\zeta({\Lambda}-1)}{{\Lambda}}}{(1+\rho\cdot\zeta{\Lambda}^{2})+\sqrt{(1-\rho\cdot\zeta{\Lambda}^{2})^{2}+4\rho{\Lambda}(\zeta-1)}}\Biggr{]}}^{1/2},\;\rho=0.35\;and\;0.05,\end{array}$ (2.8)

are bounds $c_{0.35}(\zeta,{\Lambda})<c<c_{0.05}(\zeta,{\Lambda})$ to the
exact solution of (2.1). We use these values in (2.2.1); for shorthand, we
drop the variables $(\zeta,{\Lambda})$.

  * ‘$\leqslant$’

It’s enough to have
$\tanh\big{(}{\Lambda}(X-\tanh({\Lambda}X))\big{)}\approx\zeta^{-1}c_{0.35}$,
so we set:

$c_{-}(\zeta,{\Lambda}):=\tanh\Big{(}{\Lambda}\cdot
c_{0.35}+\tanh^{-1}\big{(}\zeta^{-1}\cdot c_{0.35}\big{)}\Big{)}.$

  * ‘$\geqslant$’

It’s enough to have
$\tanh\big{(}{\Lambda}(X-\tanh({\Lambda}X))\big{)}\approx\zeta^{-1}c_{0.05}$,
so we set:

$c_{+}(\zeta,{\Lambda}):=\tanh\Biggl{(}\frac{\tanh^{-1}(\zeta^{-1}c_{0.05})+\frac{{\Lambda}\big{(}\sinh(2{\Lambda}c_{0.05})-2{\Lambda}c_{0.05}\big{)}}{2\cosh^{2}({\Lambda}c_{0.05})}}{1-\frac{{\Lambda}c_{0.05}}{\cosh^{2}({\Lambda}c_{0.05})}}\Biggr{)}.$

  * Remark 2.2

(i) One should keep in mind that $c_{+}(\zeta,{\Lambda})\leqslant c\leqslant
c_{-}(\zeta,{\Lambda})$. The notation is chosen in such a way that the $-$
(resp. $+$) indices determines lower (resp. upper) bounds for the exact
solution of the heat equation.

(ii) The equations (‘$\leqslant$’ ‣ 2.2.1), (‘$\geqslant$’ ‣ 2.2.1) involve
$c_{0}({\Lambda})$, which is implicitly defined. For numerical computations,
we replace it by $30$ iterations of $\tanh({\Lambda}\cdot)$ that is:

$c_{0}({\Lambda})\approx\underbrace{\tanh({\Lambda}\cdot\tanh(\dots(\tanh({\Lambda})\dots)))}_{30\;\text{times}},\quad{\Lambda}\in[1.5,5].$

For large values of ${\Lambda}$, one may use the estimates (2.2.2).

###  2.3. A transcendental equation II

For real numbers $z_{1}>z_{0}>0$, $L>0$, $\beta\geqslant\alpha>0$, we consider
the equation in $C\geqslant 0$,

  * $\displaystyle\frac{z_{1}-\alpha C}{z_{1}+\beta C}=\frac{z_{0}-\alpha C}{z_{0}+\beta C}\cdot e^{2LC}.$ (2.9)

(For our purposes, we’ll need $\beta=z_{0}=1$, only.) The substitutions

$\zeta_{i}:=z_{i}+\frac{\beta-\alpha}{2}C\;(i=0,1),\quad
c:=\frac{\beta+\alpha}{2\zeta_{0}}C,\quad{\Lambda}:=\frac{2\zeta_{0}L}{\beta+\alpha},\quad\zeta:=\frac{\zeta_{1}}{\zeta_{0}},$

reduce this equation to (2.1), but extra complications occur. The unknown $C$
is given by the equation

  * $\displaystyle C=\frac{2z_{0}c}{(\beta+\alpha)-(\beta-\alpha)c},$ (2.10)

so $\zeta,c$ are related. (Before, $c$ and $\zeta$ were independent.)
Fortunately, we are interested in solving the inequalities ‘$\leqslant$’ and
‘$\geqslant$’ in (2.3), which lead to

$\zeta=\frac{z_{1}+\frac{\beta-\alpha}{2}C}{z_{0}+\frac{\beta-\alpha}{2}C}=c\cdot\frac{1-c\cdot\tanh({\Lambda}c)}{c-\tanh({\Lambda}c)}.$

with inequalities in the same direction, respectively.

  * –

The left-hand side is decreasing in the variable $C$ (see (2.3)), and we have

    * $\displaystyle\frac{z_{0}c_{0}({\Lambda})}{\beta}\leqslant\frac{2z_{0}c_{0}({\Lambda})}{(\beta+\alpha)-(\beta-\alpha)c_{0}({\Lambda})}<C<\frac{z_{0}}{\alpha}.$ (2.11)

  * –

The right-hand side is increasing in ${\Lambda}$, and we have
$\frac{z_{0}}{\beta}L\leqslant{\Lambda}\leqslant\frac{z_{0}}{\alpha}L$.

Thus we obtain approximate solutions of the inequalities:

  * ‘$\leqslant$’

It requires solving

$\zeta_{-}:=\frac{z_{1}+\frac{(\beta-\alpha)z_{0}}{2\beta}c_{0}(\frac{z_{0}}{\beta}L)}{z_{0}+\frac{(\beta-\alpha)z_{0}}{2\beta}c_{0}(\frac{z_{0}}{\beta}L)}\leqslant
c\cdot\frac{1-c\cdot\tanh(\frac{z_{0}}{\beta}Lc)}{c-\tanh(\frac{z_{0}}{\beta}Lc)},$

which is satisfied by $c=c_{-}(\zeta_{-},\frac{z_{0}}{\beta}L)$.

  * ‘$\geqslant$’

It requires solving

$\zeta_{+}:=\frac{z_{1}+\frac{(\beta-\alpha)z_{0}}{2\alpha}}{z_{0}+\frac{(\beta-\alpha)z_{0}}{2\alpha}}\geqslant
c\cdot\frac{1-c\cdot\tanh(\frac{z_{0}}{\alpha}Lc)}{c-\tanh(\frac{z_{0}}{\alpha}Lc)},$

which is satisfied by $c=c_{+}(\zeta_{+},\frac{z_{0}}{\alpha}L)$.

The corresponding values $C_{\pm}$ are obtained by inserting $c_{\pm}$ into
(2.3).

###  2.4. Solution to an ODE

Some calculations yield the following:

  * Proposition 2.3

Given $0<\varepsilon<1$ and $0<C<\varepsilon^{-1}$. The solution to the ODE

$\frac{y^{\prime}}{\sqrt{2}}=\frac{2{\tilde{B}}}{1+\varepsilon}\bigl{(}y^{2.5}+(1-\varepsilon)Cy-\varepsilon
C^{2}y^{-0.5}\bigr{)},\quad y(1)=T,$

is the function

    * $\displaystyle y(x)=y_{{\tilde{B}},C,\varepsilon}(x):=\biggl{[}C\cdot\frac{T^{1.5}\cdot(1+\varepsilon\cdot e^{3\sqrt{2}{\tilde{B}}C(1-x)})+\varepsilon C(e^{3\sqrt{2}{\tilde{B}}C(1-x)}-1)}{T^{1.5}(e^{3\sqrt{2}{\tilde{B}}C(1-x)}-1)+C(e^{3\sqrt{2}{\tilde{B}}C(1-x)}+\varepsilon)}\biggr{]}^{2/3}.$ (2.12)

This function has the boundary value $y(0)=1$ precisely when
${\tilde{B}},C,\varepsilon$ satisfy the equation (2.3), with
$z_{0}=1,z_{1}=T,\alpha=\varepsilon,\beta=1$.

###  2.5. Criteria for upper/lower envelopes

We will determine global, as well as partial envelopes for $u_{\rm
heat},y_{\rm heat}$. For deciding that a function is indeed a global envelope,
we’ll use the following criterion:

  * Proposition 2.4

Let $f,g$ be continuous, positive on $[m,M]$, such that one of them is
strictly increasing. We consider the solutions to the following boundary-value
problems:

$\begin{array}[]{llr}&y_{0}^{\prime\prime}=f(y_{0}),&y_{0}(0),y_{0}(1)\;given\;in\;[m,M],\\\
&y^{\prime\prime}=g(y),&y(0),y(1)\;given\;in\;[m,M].\\\ \end{array}$

If one has the following three inequalities

$g-f,\;y_{0}(0)-y(0),\;y_{0}(1)-y(1)\geqslant 0,\quad\text{(resp. $\leqslant
0$)}$

then we have $y_{0}-y\geqslant 0$ (resp. $\leqslant 0$).

In this article, we’ll consider $f(y)=B^{2}(5y^{4}-5)$. The proposition can be
reformulated in terms of the so-called ‘residue functional’ corresponding to
it, which associates the function

${\rm resd}(y):=y^{\prime\prime}-B^{2}(5y^{4}-5),$

to any (differentiable) function $y$. The residue measures the deviation of
$y$ from being a solution to the given ODE. Our criterion reads as follows:

  * –

If $y(0)\geqslant 1,y(1)\geqslant T$, and ${\rm resd}(y)\leqslant 0$, then $y$
is greater than $y_{\rm heat}$.

  * –

If $y(0)\leqslant 1,y(1)\leqslant T$, and ${\rm resd}(y)\geqslant 0$, then $y$
is smaller than $y_{\rm heat}$.

This is a criterion for global envelopes. However, we’ll be also interested in
partial envelopes, along the interesting boundary layer.

  * Proposition 2.5

Let $F,G$ be positive, differentiable on $[1,T]$, and consider the IVPs

$\begin{array}[]{llll}y_{0}^{\prime}=\sqrt{2}\cdot
F(y_{0}),&\;y^{\prime}=\sqrt{2}\cdot
G(y),&&y_{0}(1)=y_{0}(1)=T\;\text{given}.\end{array}$

If the following inequality is satisfied (where $r\in[1,T]$):

$G-F\geqslant 0,\;\text{for}\;y\in[r,T],\quad\text{(resp. $\leqslant 0$)},$

then $y_{0}-y\geqslant 0$ (resp. $\leqslant 0$), whenever $y\in[r,T]$.

The function of interest will be $F(y)=B^{2}\cdot\sqrt{y^{5}-5y+4+\delta}$,
with $\delta\geqslant 0$. Let us explain how this will be used to determine
partial envelopes in Propositions 5.1.2 and 5.2.1. The IVPs are obtained by
respectively integrating the BVPs, so

$F(y)=\sqrt{\delta+\int_{0}^{y}f(\xi){\rm d}\xi},\;\text{ and similarly
for}\;G.$

Suppose we are in the following situation: $G(r)\geqslant F(r)$ and there are
solutions to the IVPs above, such that ${\rm resd}(y)=g(y)-f(y)\geqslant
0,\;y\in[r,T]$; this condition holds, in particular, if ${\rm
resd}(y=r)\geqslant 0$ and ${\rm resd}(y)$ in increasing on $[r,T]$). Then
Proposition 2.5 applies.

##  3\. The heat conduction equation: I

###  3.1. Estimates

We are interested in bounding the solution of the BVP:

  * $\displaystyle u^{\prime\prime}_{\rm heat}=b^{2}(u^{4}_{\rm heat}-t^{4}),\quad u_{\rm heat}(0)=1,\;u_{\rm heat}(1)=t\in(0,1).$ (3.1)

We perform several changes of variables, to bring the problem into a more
suitable form. The substitution $x\mapsto 1-x$ yields

  * $\displaystyle u^{\prime\prime}=b^{2}(u^{4}-t^{4}),\quad u(0)=t\in(0,1),\;u(1)=1.$ (3.2)

By integrating, we obtain

$\frac{(u^{\prime})^{2}}{2t^{2}}=\frac{b^{2}t^{3}}{5}\Bigl{(}\frac{u^{5}}{t^{5}}-5\frac{u}{t}+A\Bigr{)}.$

The derivative at the left-end of the interval $[0,1]$ is expected to be close
to zero. We set $A=4+\delta$ and $y:=u/t$, so the equation becomes:

  * $\displaystyle\begin{array}[]{ll}\displaystyle\frac{y^{\prime}_{\rm heat}}{\sqrt{2}}=B\cdot\sqrt{y^{5}_{\rm heat}-5y_{\rm heat}+4+\delta},&\displaystyle y_{\rm heat}(0)=1,y_{\rm heat}(1)=t^{-1},\\\\[8.61108pt] \displaystyle y^{\prime\prime}_{\rm heat}=5B^{2}(y^{4}_{\rm heat}-1),&\displaystyle B:=\frac{bt^{3/2}}{\sqrt{5}}.\end{array}$ (3.3)

Clearly, the equation can’t be integrated exactly, so we need to approximate
it. At the left-end of $[0,1]$, we have $y_{\rm heat}(0)=1$. However, at the
right-end, the function $y_{\rm heat}$ take large values, for $t$ close to
zero, so we need to take into account the large-$y$ behaviour. Therefore, we
consider the boundary value problem

${\frac{y^{\prime}}{\sqrt{2}}=B\cdot G(y),\quad y(0)=1,y(1)=\frac{1}{t},}$

with the right-hand side satisfying the following requirements:

  * –

The difference $\sqrt{y^{5}-5y+4+\delta}-G(y)$ is small, for all $y\geqslant
1$;

  * –

The equation $y^{\prime}=B\cdot G(y)$ is explicitly integrable.

For large $y$, we have the expansion

$\sqrt{y^{5}-5y+4+\delta}=y^{\frac{5}{2}}-\frac{5}{2}y^{-\frac{3}{2}}+\frac{4+\delta}{2}y^{-\frac{5}{2}}+\dots=y^{\frac{5}{2}}\cdot\biggl{(}1-\frac{5}{2}y^{-4}+\frac{4+\delta}{2}y^{-5}+\cdots\biggr{)},$

which hints to considering $G(y)=y^{5/2}\cdot c(y)$, such that the
approximation is still good for $y\approx 1$. Thus the function $c(y)$ must
satisfy $c(1)\approx 0$ and $\mathop{\lim}\limits_{y\to\infty}c(y)=1$. The
asymptotic expansion hints to considering

$G(y)=y^{2.5}\cdot(1-\text{const.}\cdot y^{-2.5-r})=y^{2.5}-\text{const.}\cdot
y^{-r},$

with $\text{const.}\approx 1$. It remains to determine a value of $r$ which is
appropriate for both requirements. The numerator of the difference

$\sqrt{y^{5}-5y+4}-(y^{5/2}-y^{-r})=\frac{-5y+4+2y^{2.5-r}-y^{-2r}}{\sqrt{y^{5}-5y+4}+(y^{5/2}-y^{-r})},$

vanishes of second order at $y=1$ for $r_{0}=\sqrt{10}-2.5\approx 0.6623$; for
approximation purposes, a convenient value for $r$ should be chosen nearby.
However, the resulting ODE must be explicitly integrable, too, and this
determines $r=0.5$.

###  3.2. A boundary value problem I

The value $r=0.5$ leads to the following BVP (where we set $T:=t^{-1}$):

  * $\displaystyle\begin{array}[]{l|l}\frac{y^{\prime}}{\sqrt{2}}=B\cdot(y^{2.5}-C^{2}y^{-0.5}),&y^{\prime\prime}=B^{2}(5y^{4}-4C^{2}y-C^{4}y^{-2}),\\\ y(0)=1,y(1)=T,&y(0)=1,y(1)=T.\end{array}$ (3.4)

In the first formulation, the constant $C>0$ is determined by the boundary
conditions. For approximating the exact solution, we change our point of view.
We drop the boundary value at $x=0$, and consider the initial value problem
(IVP):

  * $\displaystyle\begin{array}[]{l}\frac{y^{\prime}}{\sqrt{2}}={\tilde{B}}\cdot(y^{2.5}-C^{2}\cdot y^{-0.5}),\quad y\geqslant 1,\;C\in(0,1],\\\\[4.30554pt] y(1)=T,\end{array}$ (3.5)

whose general solution is:

  * $\displaystyle\begin{array}[]{l}\displaystyle\frac{T^{1.5}-C}{T^{1.5}+C}=\frac{y^{1.5}-C}{y^{1.5}+C}\cdot e^{3\sqrt{2}{\tilde{B}}C\cdot(1-x)},\quad\frac{T^{1.5}-C}{T^{1.5}+C}=\frac{1-C}{1+C}\cdot e^{3\sqrt{2}{\tilde{B}}C},\\\\[8.61108pt] \displaystyle y={\Biggl{[}\frac{C}{\tanh\Big{(}\tanh^{-1}(Ct^{1.5})+1.5\sqrt{2}{\tilde{B}}C(1-x)\Big{)}}\Biggr{]}}^{2/3}\kern-4.30554pt.\end{array}$ (3.6)

(The parameter ${\tilde{B}}$ should be thought off as an ‘improved value’ for
$B$, which increases the accuracy of the estimates.) The equation defining $C$
was analysed in Section 2, with the parameter
${\Lambda}=1.5\sqrt{2}{\tilde{B}}$; by our discussion, we have $C\leqslant 1$.
To obtain an upper bound for (3.1), the second derivative should satisfy
Proposition 2.5:

  * $\displaystyle\begin{array}[]{l}y^{\prime\prime}={\tilde{B}}^{2}(5y^{4}-4C^{2}y-C^{4}y^{-2})\leqslant B^{2}(5y^{4}-5),\quad\forall\,y\in[1,T],\quad y(0)\geqslant 1.\end{array}$ (3.7)

The left-hand side is decreasing in $C$, so it’s enough to ensure that the
inequality is satisfied for $C=1$. Let

  * $\displaystyle{\tilde{B}_{+}}:=(\tilde{q}_{+})^{-1}\cdot B,\;\text{with}\;\tilde{q}_{+}:=\Bigl{[}\frac{(1-t^{3})(5+t^{3})}{5(1-t^{4})}\Bigr{]}^{1/2}\in[\sqrt{0.9},1].$ (3.8)

  * Lemma 3.1

The function
$\displaystyle\tilde{y}_{+}(x):={\Big{[}\tanh\Bigl{(}\tanh^{-1}(t^{1.5})+1.5\sqrt{2}{\tilde{B}_{+}}(1-x)\Bigr{)}\Big{]}}^{-\frac{2}{3}},$
satisfies

$\begin{array}[]{l}\tilde{y}_{+}^{\prime\prime}={({\tilde{B}_{+}})}^{2}(5\tilde{y}_{+}^{4}-4\tilde{y}_{+}-\tilde{y}_{+}^{-2}),\\\\[8.61108pt]
\tilde{y}_{+}(1)=T=t^{-1},\;\tilde{y}_{+}(0)=\biggl{(}1+\frac{2}{\exp\bigl{(}2\tanh^{-1}(t^{1.5})+3\sqrt{2}{\tilde{B}_{+}}\bigr{)}-1}\biggr{)}^{\frac{2}{3}},\end{array}$

and it is greater than the exact solution $y_{\rm heat}$ of the BVP (3.1).

  * Proof.

The function $\tilde{y}_{+}$ satisfies the conditions of Proposition 2.5.
$\Box$

We briefly pause the analysis to analyse the behaviour of the function $y_{\rm
heat}$, in order to justify the subsequent considerations and refinements.

###  3.3. Estimating the derivative

The major issue for the shooting method [4] is the difficulty to estimate the
initial derivative $y^{\prime}_{\rm heat}(0)=\sqrt{\delta}$. A crucial by-
product of Lemma 3.2 is an upper bound for $\delta$ and, a posteriori, for the
derivative at $x=1$.

  * Proposition 3.2

The following estimates hold for the derivative of $y_{\rm heat}$:

    * $\begin{array}[]{l}\displaystyle\delta\leqslant\Delta:={\Bigg{[}\frac{\tilde{y}_{+}(0.25)-1}{0.25\sqrt{2}B}\Bigg{]}}^{2}\leqslant\frac{15}{B^{2}}\cdot\frac{1}{(e^{3B}-1)^{2}}\Big{(}1+\frac{1}{e^{3B}-1}\Big{)}^{2};\\\ \displaystyle 0\leqslant y^{\prime}_{\rm heat}(1)-\sqrt{2}B\sqrt{T^{5}-5T+4}\leqslant\frac{\sqrt{2}B\delta}{2\sqrt{T^{5}-5T+4}}\leqslant\frac{11\Big{(}1+\frac{1}{e^{3B}-1}\Big{)}^{2}}{B(e^{3B}-1)^{2}\cdot\sqrt{T^{5}-5T+4}}.\end{array}$ (3.9)

The moral is that, already for $B\geqslant 2$, the derivative $y_{\rm
heat}^{\prime}(1)$ differs from the (universal) lower bound
$\sqrt{2}B\sqrt{T^{5}-5T+4}$ by a tiny amount, of size
$o(B^{-1}e^{-6B}T^{-2.5})$. Indeed, for $B\geqslant\frac{1}{3}$, one has
$\delta\leqslant 100B^{-2}e^{-6B}.$

  * Proof.

Note that $y_{\rm heat}$ is convex —that is, $y_{\rm heat}^{\prime\prime}>0$—,
so its derivative at $x=0$ is less than the slope of the secant line between
$x=0$ and $x=h=0.25$. Since $y_{\rm heat}$ is less than $\tilde{y}_{+}$, we
deduce the following estimate:

$\begin{array}[]{rl}\sqrt{2}B\sqrt{\delta}=y_{\rm
heat}^{\prime}(0)&\leqslant\frac{\tilde{y}_{+}(h)-1}{h}\leqslant\frac{\tilde{y}_{+}(h)^{3}-1}{3h}=\frac{2\big{(}\tilde{y}_{+}(h)^{3/2}-1\big{)}+\big{(}\tilde{y}_{+}(h)^{3/2}-1\big{)}^{2}}{3h}.\end{array}$

By inserting $x=h$ into the defining formula, we obtain
$\tilde{y}_{+}(h)^{3/2}-1\leqslant\frac{2}{e^{3.1B}-1}$, thus

$\sqrt{\delta}\leqslant\frac{4}{h\cdot 3\sqrt{2}\cdot
B}\cdot\frac{1}{e^{3.1B}-1}\Big{(}1+\frac{1}{e^{3B}-1}\Big{)}.$

It remains to use this estimate, to bound the derivative of $y_{\rm heat}$ at
$x=1$. $\Box$

###  3.4. Application to the ‘shooting’ method

The content of this section should probably belong to Section 7, which is
devoted to numerical simulations. However, we believe that this short
digression is useful to present some applications following from theoretical
considerations.

The shooting method is doomed to fail without an accurate knowledge of the
initial derivative of $y_{\rm heat}$. According to Proposition 3.3, we have

$0\leqslant\delta\leqslant\Delta={\Bigg{[}\frac{\tilde{y}_{+}(0.25)-1}{0.25\sqrt{2}B}\Bigg{]}}^{2}.$

This bare inequality already yields three possible ways of approximating
$y_{\rm heat}$ (the abbreviation ‘RK’ stands for Runge-Kutta):

  1. (i)

Let $\delta=0$, and consider the IVP
$y^{\prime}/\sqrt{2}=B\sqrt{y^{5}-5y+4},\;y(1)=T$, whose solution is denoted
by $y_{RK+}$. Its second derivative satisfies the same equation
$y^{\prime\prime}=5B^{2}(y^{4}-1)$ as $y_{\rm heat}$, and
$y_{RK+}^{\prime}(1)<y_{\rm heat}^{\prime}(1)$, hence $y_{RK+}>y_{\rm heat}$.

  2. (ii)

Let $\delta=\Delta$, and consider
$y^{\prime}/\sqrt{2}=B\sqrt{y^{5}-5y+4+\Delta},\;y(1)=T$, whose solution is
denoted by $y_{RK-}$. Its second derivative satisfies
$y^{\prime\prime}=5B^{2}(y^{4}-1)$, same as $y_{\rm heat}$, and
$y_{RK-}^{\prime}(1)>y_{\rm heat}^{\prime}(1)$, hence $y_{RK-}<y_{\rm heat}$.

  3. (iii)

Let
$\delta=d:={\bigg{[}\frac{\tilde{y}_{+}(0.25)-\tilde{y}_{+}(0)}{0.25\sqrt{2}B}\bigg{]}}^{2}$,
which is between $0$ and $\Delta$, and consider the corresponding IVP. Its
solution $y_{RK}$ will approximate $y_{\rm heat}$ the best.

One is led to asking what is the precision of these approximations.

  * Corollary 3.3

The error of the approximation $y_{\rm heat}(x)\approx y_{RK}(x)$ is at most
$err_{y}(x)=y_{RK+}(x)-y_{RK-}(x)$. The size of the maximal error equals
$O(e^{-3B})$.

The error for $u_{\rm heat}(1-x)\approx t\cdot y_{RK}(x)$ has the order
$O(t\cdot e^{-3B})$.

  * Proof.

The error function is decreasing, has negative derivative, so the maximum is
attained at $x=0$, and is bounded above by $const\cdot\sqrt{\Delta}$. $\Box$

We probe this matter numerically, by plotting a few examples with the Runge-
Kutta method (using MAPLE).

Table 1. Examples of plots using the Runge-Kutta method values of parameters | graphs of envelopes | graph of error
---|---|---
|  $y_{RK+}\,$(dots), $y_{RK-}\,$(line) | $err_{y}(x)=y_{RK+}-y_{RK-}$
$\begin{array}[]{l}b=10,\\\ t=0.3,\\\ B=0.73,\\\ max.err_{y}=0.035,\\\ \frac{\ln(max.err_{y})}{B}=-4.55.\end{array}$ |  |
$\begin{array}[]{l}b=55,\\\ t=0.1,\\\ B=0.78,\\\ max.err_{y}=0.038,\\\ \frac{\ln(max.err_{y})}{B}=-4.18.\end{array}$ |  |
$\begin{array}[]{l}b=30,\\\ t=0.7,\\\ B=7.85,\\\ max.err_{y}=4.7\cdot 10^{-11},\\\ \frac{\ln(max.err_{y})}{B}=-3.02.\end{array}$ |  |

The conclusion is that the Runge-Kutta method can be successfully run, as long
as there is enough computational power available; high values of the
parameters require long computational time. It’s certainly remarkable that the
error decreases exponentially fast with $B$. (The last graphs of the envelopes
overlap.)

###  3.5. The boundary layer condition

The estimates (3.3) are essential for understanding the behaviour of the
solution to the equation (3.1), especially the boundary layer phenomenon.

Note that the secant line between $(0,1)$ and $(1,T)$ is above the graph of $y_{\rm heat}$, so $y_{\rm heat}^{\prime}(1)>T-1$. Now consider the angles, denoted $a$ and $s$ in the figure. For a parameter $\kappa>1$, the $\kappa$-boundary layer property should be that $\displaystyle a<\frac{s}{\sqrt{2}\kappa}=\frac{\tan^{-1}\bigl{(}\frac{1}{T-1}\bigr{)}}{\sqrt{2}\kappa}.$ (3.10) |
---|---

For $T$ —same for $t=T^{-1}$— close to $1$, the function $y_{\rm heat}$ is
approximately constant $1$. Thus, we are interested in larger values, such as
$T\geqslant 1.5$ that is $t\leqslant 2/3$. Since $y_{\rm
heat}^{\prime}(1)=a^{-1}$ is approximately $\sqrt{2}BT^{2.5}$, we replace the
inequality above with $y_{\rm heat}^{\prime}(1)>\sqrt{2}\kappa T$, which
yields:

  * $\displaystyle BT^{1.5}>\kappa\quad\Leftrightarrow\quad b>\sqrt{5}\kappa.$ (3.11)

We choose the value $\kappa=35$ for our parameter: then (3.5) implies, for
$T>2$, that the angle $a$ is less than $1^{\circ}$. From a visual standpoint,
at least for $x$ very close to $1$, this should be steep enough to claim the
existence of a “boundary layer”.

Now, let us compare this geometric approach with the analytic data: we plot
the graph of $\tilde{y}_{+}$ for several values of the parameters $B$ and $T$.
Since $\tilde{y}_{+}$ is greater than $y_{\rm heat}$, the graph of latter is
below the former, thus $y_{\rm heat}$ is “at least as horizontal” about $x=0$
and “at least as vertical” about $x=1$. We verify whether, from a visual
viewpoint, the boundary layer property is imposed only by the size of the
derivative at $x=1$. The will see that we must slightly adjust our previous
considerations.

The estimate for $y_{\rm heat}^{\prime}(0)$ shows that already for low values of $B$ (e.g. $B\geqslant 0.5$), the graph of $y_{\rm heat}$ starts off almost horizontally at $x=0$. This is not true for low values of $B$! |
---|---

On the right, $BT^{1.5}>35$ is large, yet there is no (visual) boundary layer. (The plot is only over the interval $x\in[0.9,1]$). The explanation for this matter is the following: the condition $y_{\rm heat}^{\prime}(1)\gg 0$ doesn’t necessarily imply that the graph “falls vertically” along the line $x=1$. The idea which emerges is that, to ensure steep-falling, one should impose that $y_{\rm heat}^{\prime}$ is large at a suitable value $x<1$. (Note that, since $y_{\rm heat}^{\prime}$ is increasing, $y_{\rm heat}^{\prime}(1)$ is large, too.) |
---|---

  * Definition 3.4

We say that $y_{\rm heat}$ possesses the boundary layer property if the
following condition holds:

$B\cdot T^{1.25}\geqslant 50\quad\text{or}\quad b\sqrt[4]{t}>110.$

  * –

Let $h:=y_{\rm heat}^{-1}(\sqrt{T})<1$. The estimate (3.3) shows that our
definition basically means

$y_{\rm heat}^{\prime}(h)\geqslant 50\sqrt{2}>\tan(89^{\circ})\approx 57.29,$

so the angle $a$ in (3.5) is less than $1^{\circ}$, too.

  * –

Our benchmark for declaring the “boundary layer” was the angle $1^{\circ}$.
One obtains steeper boundary layer —the graph of $y_{\rm heat}$ is even more
vertical— by increasing the lower bound of $BT^{1.25}$.

We conclude with another geometric description of the boundary-layer property.

  * Corollary 3.5

The existence of the boundary layer implies that the function dramatically
decreases over a very short interval. The ratio between the variation of
$y_{\rm heat}$ over the interval $[1-\sqrt{5}/b,1]$ and its total variation
over $[0,1]$ is at least $0.5$.

  * Proof.

We observe that $\xi=\sqrt{5}/b=t^{1.5}/B\leqslant t^{0.25}/50\leqslant 0.02$,
so the indicated interval is indeed short. Moreover, since
$\tilde{y}_{+}\geqslant y_{\rm heat}$ and ${\tilde{B}_{+}}\geqslant B$, we
have:

$\frac{T-y_{\rm
heat}(1-\xi)}{T-1}\geqslant\frac{T-\tilde{y}_{+}(1-\xi)}{T-1}\geqslant\frac{1-t\,\cdot\,{\tanh(\,\tanh^{-1}(t^{1.5})+1.5\sqrt{2}t^{1.5})\,)^{-\frac{2}{3}}}}{1-t}>0.5.$

The last inequality is obtained by plotting the function (for $t\in[0,1]$) and
reading off its minimum ($\approx 0.53$). $\Box$

###  3.6. Inconveniences

We justify the necessity for further investigation.

  * –

First, one can’t hope to obtain the opposite inequality (3.2) —it’s clear for
$T\gg 0$—, so the IVP (3.2) is not suitable for deducing a lower bound to
$y_{\rm heat}$.

  * –

Second, the sensitive dependence of the problem on the parameters $b,t$, makes
that $\tilde{y}_{+}$ above is typically a rather loose estimate for $y_{\rm
heat}$.
In Table LABEL:tab:y+res, we plotted $\tilde{y}_{+}$ and its residue

${\rm
resd}(\tilde{y}_{+}):={{\tilde{B}_{+}}}^{2}(5\tilde{y}_{+}^{4}-4\tilde{y}_{+}-\tilde{y}_{+}^{-2})-B^{2}(5\tilde{y}_{+}^{4}-5),$

for $B=13,T=3$. We are in the boundary layer case $BT^{1.25}\approx 51$, the
graph falls vertically, but the residue is large near $x=1$; its $L^{2}$-norm
is $\|{\rm resd}(\tilde{y}_{+})\|_{L^{2}{([0,1]})}\approx 36.7$.

Table 2. Graphs of $\tilde{y}_{+}$ and ${\rm resd}(\tilde{y}_{+})$. |
---|---

##  4\. The heat conduction equation: II

Now we refine the previous methods.

###  4.1. A boundary value problem II

The change of variables $z=y^{1.5}$ in (3.2) yields the ODE
$2z^{\prime}=3\sqrt{2}B(z^{2}-C^{2})$, which requires integrating

$\frac{2{\rm d}z}{z^{2}-C^{2}}=\frac{{\rm
d}z}{C}\Big{[}\frac{1}{z-C}-\frac{1}{z+C}\Big{]}.$

Here we consider a one-parameter deformation of this $1$-form, namely:

$\frac{{\rm
d}z}{C}\Big{[}\frac{1}{z-C\varepsilon}-\frac{1}{z+C}\Big{]}=\frac{(1+\varepsilon){\rm
d}z}{z^{2}+(1-\varepsilon)Cz-\varepsilon C^{2}},$

where $\varepsilon<1$ is introduced to shift the pole of the initial $1$-form
away from $z=C\approx 1$, since poles yield major computational errors. The
function $z$ solves

$\frac{{\rm
d}z}{C}\Big{[}\frac{1}{z-C\varepsilon}-\frac{1}{z+C}\Big{]}=3\sqrt{2}B{\rm
d}x,\quad z(0)=z_{0}\approx 1,\;z(1)=T^{1.5},$

so it satisfies the equation

  * $\displaystyle\frac{z-\varepsilon C}{z+C}=\frac{T^{1.5}-\varepsilon C}{T^{1.5}+C}\cdot e^{-3\sqrt{2}BC(1-x)},$ (4.1)

with $C$ determined by the equality:

  * $\displaystyle\frac{T^{1.5}-\varepsilon C}{T^{1.5}+C}=\frac{z_{0}-\varepsilon C}{z_{0}+C}\cdot e^{2LC},\quad L:=1.5\sqrt{2}B.$ (4.2)

We recover the situation described in Remark 2.2.3.

For the initial function $y=z^{2/3}$, we are considering the solution to the
new BVP,

  * $\frac{y^{\prime}}{\sqrt{2}}=\frac{2{\tilde{B}}}{1+\varepsilon}\bigl{(}y^{2.5}+(1-\varepsilon)Cy-\varepsilon C^{2}y^{-0.5}\bigr{)},\quad y(0)=1,\;y(1)=T,$ (4.3)

whose second derivative satisfies:

  * $y^{\prime\prime}=\frac{4{\tilde{B}}^{2}}{(1+\varepsilon)^{2}}\Bigl{(}5y^{4}+7(1-\varepsilon)Cy^{2.5}-2(4\varepsilon-1-\varepsilon^{2})C^{2}y-\varepsilon(1-\varepsilon)C^{3}y^{-0.5}-\varepsilon^{2}C^{4}y^{-2}\Bigr{)}.$ (4.4)

The discussion in Section 3.2 corresponds to $\varepsilon=C=1$. We introduce
the notation:

$B=q\cdot{\tilde{B}}.$

The residue function of $y$ satisfies the formula:

$\begin{array}[]{rl}\frac{{\rm
resd}(y)}{B^{2}}&=R(y,q,C,\varepsilon)=\bigl{(}1-q^{2}\frac{(1+\varepsilon)^{2}}{4}\bigr{)}\cdot
5y^{4}+7(1-\varepsilon)Cy^{2.5}-2(4\varepsilon-1-\varepsilon^{2})C^{2}y\\\\[8.61108pt]
&-\varepsilon(1-\varepsilon)C^{3}y^{-0.5}-\varepsilon^{2}C^{4}y^{-2}+5q^{2}\frac{(1+\varepsilon)^{2}}{4},\end{array}$

involving the parameters $q,\varepsilon,C$, which must be finely tuned so that
Proposition 2.5 will be eventually satisfied: the solution to the BVP (4.1) is
a lower (resp. upper) bound for the solution of (3.1) if
$R(y,q,C,\varepsilon)$ is positive (resp. negative) for $y\in[0,1]$.

  * Procedure 4.1

    1. (i)

The value of $q$ is determined by the derivative condition

$\frac{y^{\prime}(1)}{\sqrt{2}}\approx B\sqrt{T^{5}-4T+(4+\delta)},$

and by ensuring the correct inequality; the bounds (3.3) for $\delta$ are
essential.

       * –

for $\leqslant$, we impose $y^{\prime}(1)/\sqrt{2}\leqslant
q{\tilde{B}}\cdot\sqrt{T^{5}-4T+4}$.

       * –

for $\geqslant$, we impose $y^{\prime}(y)/\sqrt{2}\geqslant
q{\tilde{B}}\cdot\sqrt{y^{5}-4y+(4+\Delta)}$, for $y=T,\sqrt{T}$;

    2. (ii)

The deformation parameter $\varepsilon$ will be typically close to $1$, and
it’s determined by the curvature/residue condition $R(y,\dots)=rB^{2},$ where
$r$ is a small real number (positive, resp. negative).

    3. (iii)

Finally, $C$ is determined by the boundary condition $y(0)\approx 1$ (cf.
Section 2); one expects it to be close to $1$, too. Note that $C$ is an
explicit function of $q,\varepsilon$.

###  4.2. Analysis of the residue

The sign of the residue function is a decisive matter. In explicit situations,
the simplest way to settle the issue is by plotting the function. Here we
analyse the sign of this expression (dependence on parameters) by relying on
analytical tools. The polynomial in four variables above has the following
properties:

  * –

$R(1,1,1,1)=0,\;\frac{{\rm d}R}{{\rm d}C}(1,1,1,1)<0,\;\frac{{\rm d}R}{{\rm
d}\varepsilon}(1,1,1,1)<0;$

  * –

$R(1,1,C,\varepsilon)=C^{2}(2+C-C^{2})\varepsilon^{2}-C(7+8C+C^{2})\varepsilon+5+7C+2C^{2};$

  * –

The quadratic equation $R(1,1,1,\varepsilon)=0$ admits two roots: $1$ and $7$.
For small variations of the parameters, these roots will persist, and we
denote by $\varepsilon(y,q,C,r)$ the root of $R(y,q,C,\varepsilon)=r$ which is
close to $1$. For shorthand, we abusively denote it by $\varepsilon(y)$, since
we are mainly interested in the dependence in $y$.

  * Lemma 4.2

We denote $s:=y^{-1}\in[0,1]$, and let:

$\tilde{\varepsilon}_{-}(s):=1-\frac{4}{5}s^{3}+\frac{3}{5}s^{4},\qquad\tilde{\varepsilon}_{+}(s):=1-\frac{4}{5}s^{3}+\frac{8-3s}{5}s^{4}.$

Then the following statements hold:

    * –

$R(s^{-1},1,1,\tilde{\varepsilon}_{-}(s)),R(s^{-1},1,1,\tilde{\varepsilon}_{+}(s))$
are continuous, for $s\in[0,1]$.

    * –

They satisfy the inequalities:

$\forall\,s\in[0,1],\qquad R(s^{-1},1,1,\tilde{\varepsilon}_{-}(s))\geqslant
2,\quad R(s^{-1},1,1,\tilde{\varepsilon}_{+}(s))\leqslant-\frac{2}{3}.$

Thus, for small values of the parameters $q,C,r$, the root $\varepsilon(s)$ of
$R(s^{-1},q,C,\varepsilon)=r$ belongs to the interval
$[\tilde{\varepsilon}_{-}(s),\tilde{\varepsilon}_{+}(s)]$.

One should still clarify what’s the meaning of ‘small’ in the statement. The
series expansion of the root of $R(s^{-1},1,C,r)=0$ is
$\varepsilon(y)=1-\frac{4}{5}C^{2}s^{3}+\Big{(}1-\frac{r}{5}\Big{)}s^{4}+O(s^{5})$.
One may check that the Lemma still holds for $1-C\leqslant t^{3}/10$ and
$|r|\leqslant 0.5$.

  * Proof.

The continuity reduces to the fact that the series development of $R$ contains
only positive powers of $s$. The indicated inequalities are obtained by
plotting the corresponding graphs.

For the second statement, a computation shows that the power series expansion
of the root of $R(s^{-1},1,C,r)=0$ is
$\varepsilon(y)=1-\frac{4}{5}C^{2}s^{3}+\Big{(}1-\frac{r}{5}\Big{)}s^{4}+O(s^{5})$.
$\Box$

Lemma 4.2 is important to understand the monotonicity of the residue function.
Below we compute the components of the gradient of $R$; the abbreviation
${l.o.t.}$ stands for ‘lower order terms’. Also, we think off $C$ as function
of $q$ and $\varepsilon$, since $C$ is determined by the boundary condition at
$x=0$ (see Procedure 4.1); a very rough approximation would be
$\tanh(1.5\sqrt{2}Bq^{-1})$.

$\begin{array}[]{ll}\displaystyle\frac{{\rm d}R}{{\rm
d}y}&=\Big{(}1-q^{2}\frac{(1+\varepsilon)^{2}}{4}\Big{)}\cdot
20y^{3}+17.5(1-\varepsilon)Cy^{1.5}-2(4\varepsilon-1-\varepsilon^{2})C^{2}+{l.o.t.};\\\\[8.61108pt]
\displaystyle{\Bigl{.}\frac{{\rm d}R}{{\rm
d}q}\Bigr{|}}_{y=T}&=-10q\frac{(1+ee)^{2}}{4}T^{4}+7(1-\varepsilon)\frac{{\rm
d}C}{{\rm
d}q}T^{2.5}+{l.o.t.}\approx-10T^{4}-14(1-\varepsilon)BT^{2.5};\\\\[8.61108pt]
\displaystyle{\Bigl{.}\frac{{\rm d}R}{{\rm
d}\varepsilon}\Bigr{|}}_{y=T}&=-5q^{2}\frac{1+\varepsilon}{2}T^{4}-7CT^{2.5}+7(1-\varepsilon)\frac{{\rm
d}C}{{\rm d}\varepsilon}T^{2.5}+{l.o.t.}\approx-5T^{4}.\end{array}$

For $1-\varepsilon\approx\frac{4}{5}t^{3}$, as suggested by the Lemma,
$\frac{{\rm d}R}{{\rm d}y}$ is mostly positive, so $y\mapsto{\rm resd}(y)$
will be increasing almost over the whole $[0,T]$. (Here we think off $y$ as a
variable, rather than a function of $x$.) Thus, for $y=y(x)$ an approximate
solution of the heat equation, ${\rm resd}(y(x))$ increases dramatically, for
$x$ near $1$.

Furthermore, the dominant component of the gradient is $\frac{{\rm d}R}{{\rm
d}q}$; the factor of $T^{4}$ is responsible for large variations of the
residue for tiny variations of the parameters. The only way to ascertain the
increase/decrease of $R$ is by changing $q,\varepsilon$ along the gradient
flow.

##  5\. Bounding the exact solution $y_{\rm heat}$

###  5.1. A minorant $y_{-}$ of $y_{\rm heat}$

(i) determine $q$ In the sequel, we let

$M(y^{-1}):=\sqrt{1-5y^{-4}+(4+\Delta)y^{-5}},\;q(y,\varepsilon,C):=\frac{2(1+(1-\varepsilon)C\cdot
y^{-1.5}-\varepsilon C^{2}\cdot y^{-3})}{(1+\varepsilon)\cdot M(y^{-1})}.$

The slope inequality condition $y_{\rm
heat}(y)\leqslant\frac{y^{\prime}(y)}{\sqrt{2}}$ is satisfied as soon as
$q\leqslant q(y,\varepsilon,C)$; the right-hand side is decreasing in $C$, so
we may replace it by its upper bound $\varepsilon^{-1}$ (see equation (– ‣
2.3)). Hence the slope inequality at $y$ is satisfied for

$q\leqslant q(y,\varepsilon,\varepsilon^{-1}).$

####  5.1.1. A global minorant

For $\varepsilon=1,q=1$, the residue is negative; for $\varepsilon=0,q=1$, the
residue is positive. We need values of $\varepsilon,q$ in between, for which
the residue positive, for all $y\in[0,1]$, see Proposition 2.5.

  * Proposition 5.1

Let $b>0,t\in(0,1)$, $\tilde{\varepsilon}_{-}:=0.75$. Define

$\tilde{q}_{-}:=\min\Bigl{(}1.1,q(T,\tilde{\varepsilon}_{-},\tilde{\varepsilon}_{-}^{-1})\Bigr{)},\;{\tilde{B}}_{-}:=\tilde{q}_{-}^{-1}B.$

We define $\tilde{C}_{-}:=C_{-}(\tilde{q}_{-},\varepsilon_{-})$ as in Section
2. Then the function
$\tilde{y}_{-}:=y_{{\tilde{B}}_{-},C_{-},\varepsilon_{-}}$ defined by equation
(2.4) is a global lower envelope for $y_{\rm heat}$.

  * Proof.

We prove that the residue function is positive that is,
$R(y,q_{-},C,\varepsilon_{-})>0.$ Since $R$ is decreasing in $C$, we may
replace it by an upper bound; by (2.3), we have $C\leqslant\varepsilon^{-1}$.
By the previous discussion, $q$ must satisfy the inequality:

$q\leqslant
q(T,\varepsilon_{-},\varepsilon_{-}^{-1})\leqslant\frac{2(1+(1-\varepsilon_{-})C\cdot
t^{1.5}-\varepsilon_{-}C^{2}\cdot t^{3})}{(1+\varepsilon_{-})\cdot M(t)}.$

Since $\tilde{q}_{-}\leqslant 1.1$, too, the residue condition becomes
$R(y,1.1,\tilde{\varepsilon}_{-}^{-1},\tilde{\varepsilon}_{-})\geqslant 0$,
for $y\in[1,T]$. The only negative trouble-making term is
$2(4\tilde{\varepsilon}_{-}-1-\tilde{\varepsilon}_{-}^{2})\tilde{\varepsilon}_{-}^{-1}y=\frac{23}{6}y$,
and it’s dominated either by the first two terms (for $y>1.2$) or by the last
one (for $y\in[1,1.2]$). Alternatively, one may test the inequality
graphically. $\Box$

####  5.1.2. A partial minorant

When attempting to improve the precision of the lower bound, one inevitably
faces the strong increasing tendency of $y\mapsto R(y,\dots)$, as $y$
approaches $T$, which practically makes impossible to have a global lower
bound with small residue.

The way out is to determine a partial lower bound along the boundary layer
$\sqrt{T}\leqslant y\leqslant T$, where is concentrated the interesting
information. We implement this strategy by imposing positive curvature and
slope conditions at $y=\sqrt{T}$. They yield the values of $\varepsilon$ and
$q$, respectively; since $R$ is increasing in $y$, the same conditions will
hold along $[\sqrt{T},T]$. Concretely, we proceed as follows:

  * –

Our relevant functions $q,R$ are decreasing in $\varepsilon,C$, so we fix the
‘worst’ maximal possible values:
$\varepsilon_{0}=\tilde{\varepsilon}_{+}(t),\;C_{0}:={\tilde{\varepsilon}_{-}(t)}^{-1}$.
They are needed for writing the slope and curvature conditions at
$y=\sqrt{T}$:

$\varepsilon_{1}:=\varepsilon_{0},\;q_{1}:=q(\sqrt{T},\varepsilon_{0},C_{0}),\;C_{1}:=C_{-}(\varepsilon_{1},q_{1}).$

Recall from Lemma 4.2 that $R$ is negative for
$\varepsilon=\varepsilon_{1},y=T$ —thus also for $y=\sqrt{T}$—, so the true
value, making $R$ positive, is very slightly less than that. Therefore the
slope condition, imposed by the value of $q_{1}$, will be satisfied by the
true value, too.

  * –

The precision increases if $(\varepsilon_{1},q_{1})$ is modified in the
direction of the gradient of $R$.

  * Proposition 5.2

Let $(\varepsilon_{-},q_{-},C_{-}):=(\varepsilon_{1},q_{1},C_{1})$, and
$B_{-}:=q_{-}^{-1}B$. The function

$y_{-}(x):=\biggl{[}C_{-}\cdot\frac{T^{1.5}\cdot(1+\varepsilon_{-}e^{3\sqrt{2}B_{-}C_{-}(1-x)})+\varepsilon_{-}C_{-}(e^{3\sqrt{2}B_{-}C_{-}(1-x)}-1)}{T^{1.5}(e^{3\sqrt{2}B_{-}C_{-}(1-x)}-1)+C_{-}(e^{3\sqrt{2}B_{-}C_{-}(1-x)}+\varepsilon_{-})}\biggr{]}^{2/3}$

satisfies the differential equation (4.1), and is a partial lower envelope of
the exact solution $y_{\rm heat}$, for values $y\in[\sqrt{T},T]$.

  * Proof.

We apply Proposition 2.5 and the discussion following it. $\Box$

###  5.2. An majorant $y_{+}$ for $y_{\rm heat}$

A global upper envelope was already constructed in Proposition 3.3.

####  5.2.1. A partial majorant

We construct a partial upper envelope along the boundary layer. The
monotonicity of $y\to R(y,\dots)$ implies that the residue function is
automatically negative about $x=1$, as soon as ${\rm resd}(y=T)\leqslant 0$.

We repeat the previous steps, but with the opposite inequalities; this task
will be more delicate than before. Recall that, for $q=C=\varepsilon=1$, the
function $\tilde{y}_{+}$ at Section 3.2 is an upper bound for $y_{\rm heat}$;
its major shortcoming was the (negative) magnitude of the residue about $x=1$.
To fix this issue, we need more accurate derivative and curvature conditions
at $x=1$, by following the gradient of $R$. To implement this desideratum, we
need a starting values for $q,\varepsilon,C$. The function $\tilde{y}_{+}$ is
defined (cf. equation (3.2)) using the factor

$\tilde{q}_{+}=\Bigl{[}\frac{(1-t^{3})(5+t^{3})}{5(1-t^{4})}\Bigr{]}^{1/2}=1-\frac{2}{5}t^{3}+\frac{1}{2}t^{4}-\dots$

For this reason, we take

$q_{0+}:=1-\frac{3}{5}t^{3}$

as the ‘worst’, smallest possible $q$-factor, which must be increased. For the
‘worst’ $\varepsilon$-value, Lemma 4.2 readily implies that

$\varepsilon_{0+}:=1-\frac{4}{5}t^{3}$

is quite an accurate choice. Finally, we need an (absolute) lower bound for
$C$, and this is $C_{0+}:=c_{0}(L)$. Note that $c_{0}(L)$ is defined as an
infinite iteration of the $\tanh$-function; for numerical calculations, we’ll
define $C_{0+}$ as $30$ iterations of $\tanh(0.99L)$:

$C_{0+}:=\tanh(0.99L\cdot\tanh(0.99L\cdot(\dots\cdot\tanh(0.99L)\dots))).$

(i) determine $\varepsilon$ We let

$\varepsilon_{+}:=\varepsilon(T,q_{0+},C_{0+},0)$

that is, $\varepsilon_{+}$ satisfies $R(T,q_{0+},C_{0+},\varepsilon_{+})=0$.
Since $q_{0+}$ is the smallest possible $q$-factor and $R$ is decreasing in
$q$, we’ll have $R(T,q,C_{0+},\varepsilon_{+})<0$ for our subsequently
computed values of $q$.

(ii) determine $q$ The derivative condition at $x=1$ is
$\frac{y^{\prime}(1)}{\sqrt{2}}\leqslant B\sqrt{T^{5}-5T+4}<\frac{y_{\rm
heat}^{\prime}(1)}{\sqrt{2}},$ so we must have

$q\geqslant
q_{+}(\varepsilon,C):=\frac{2}{1+\varepsilon}\cdot\frac{T^{2.5}+(1-\varepsilon)CT-\varepsilon
C^{2}T^{-0.5}}{\sqrt{T^{5}-5T+4}}.$

Since the right-hand side is decreasing with $\varepsilon,C$, we define

$q_{1+}:=q_{+}(\varepsilon_{0+},C_{0+}).$

This choice ensures that the approximate function is greater than $y_{\rm
heat}$, about $x=1$. We remark that it’s not possible to neglect
$C,\varepsilon$; doing so leads to severely wrong calculations.

(iii) determine $C$ The last condition we impose is $y_{+}(0)\geqslant 1$. The
formula (4.1) shows that this boils down to: $\frac{T^{1.5}-\varepsilon
C}{T^{1.5}+C}\geqslant\frac{1-\varepsilon C}{1+C}\cdot e^{2LC}.$ Our
discussion in Section 2.3 yields:

  * $\displaystyle\begin{array}[]{l}L_{1+}:=q_{1+}^{-1}L,\quad\zeta_{1+}:=\frac{T^{1.5}+\frac{1-\varepsilon_{1+}}{2\varepsilon_{1+}}}{1+\frac{1-\varepsilon_{1+}}{2\varepsilon_{1+}}},\\\\[8.61108pt] c_{1+}:=c_{0}(\zeta_{1+},\varepsilon_{1+}^{-1}L_{1+}),\quad C_{1+}:=C_{+}(\varepsilon_{1+},c_{1+}).\end{array}$ (5.1)

(iv) iterations The estimates greatly improve by iterating the process. The
value $\varepsilon_{+}$ above is quite precise, so we keep it. One improves
$q$ by using the gradient of $R$,

$q_{2+}:=q_{1+}+\frac{R(T,q_{1+},C_{1+},\varepsilon_{+})}{20\cdot
B^{2}}t^{4},$

and then $C_{2+}$ is determined by $\varepsilon_{+},q_{2+}$, as in (5.2.1)
above. The strange-looking modification for $q$ is motivated by the estimates
for $\frac{{\rm d}R}{{\rm d}q}$.

  * Proposition 5.3

Let $\varepsilon_{+}$ be as above and denote by $q_{+}$ and $C_{+}$ the values
obtained by two iterations; let $B_{+}:=q_{+}^{-1}B$. The function

$y_{+}(x):=\biggl{[}C_{+}\cdot\frac{T^{1.5}\cdot(1+\varepsilon_{+}e^{3\sqrt{2}B_{+}C_{+}(1-x)})+\varepsilon_{+}C_{+}(e^{3\sqrt{2}B_{+}C_{+}(1-x)}-1)}{T^{1.5}(e^{3\sqrt{2}B_{+}C_{+}(1-x)}-1)+C_{+}(e^{3\sqrt{2}B_{+}C_{+}(1-x)}+\varepsilon_{+})}\biggr{]}^{2/3}$

is a partial upper bound of the exact solution $y_{heat}$, for values
$y\in[\sqrt{T},T]$.

###### Proof.

We apply Proposition 2.5 and the discussion following it. ∎

##  6\. Back to the heat conduction equation

This study is motivated by the heat conduction equation (3.1). The bounds for
$y_{\rm heat}$ immediately yield bounds for $u_{\rm heat}$: we just set
$u=t\cdot y(1-x)$ everywhere. Thus global/partial upper/lower envelopes for
$y_{\rm heat}$ determine, respectively, the same type of envelopes for $u_{\rm
heat}$. The error for estimating $u_{\rm heat}$ is $t$ times the error for
$y_{\rm heat}$, thus smaller. The maximal computational errors for $y_{\rm
heat},u_{\rm heat}$ are compared in Table 7, for various values of $b,t$. The
Runge-Kutta method (see Section 3.4) has precision $O(t\cdot e^{-3B})$.

###  6.1. Case t = 0

An issue which naturally raises is to determine an approximate solution of the
initial heat equation, for the limiting value $t=0$. In physical sense, this
situation occurs when the ambient (final) temperature is negligible compared
to the temperature of the source. Note that $B=O(t^{1.5})$, it becomes small
for $t$ approaching $0$, so this situation is exactly opposite to the ‘$B$
large’ case.

  * Proposition 6.1

An upper envelope for the BVP,

$\;u_{0}^{\prime\prime}=b^{2}u_{0}^{4},\;u_{0}(0)=1,\;u_{0}(1)=0,\;$

is

$\displaystyle\tilde{u}_{0+}(x):={\Big{[}1+\frac{1.5\sqrt{2}}{\sqrt{5}}bx\Big{]}}^{-2/3}.$

The largest deviation from the exact solution of the BVP —the maximal error—
is attained at $x=0$, and decreases with $b$:

$max.err_{\tilde{u}_{0+}}={\Big{[}1+\frac{1.5\sqrt{2}}{\sqrt{5}}b\Big{]}}^{-2/3}\approx{\Big{[}1+0.95\cdot
b\Big{]}}^{-2/3}.$

Numerical data in the next section shows that the global upper bound
${\tilde{u}}_{+}$ for $y_{\rm heat}$ is quite precise already. Below, we’ll
see a convenient method to determine accurate lower envelope, too.

  * Proof.

Just compute $\displaystyle\lim_{t\to 0^{+}}t\cdot\tilde{y}_{+}(1-x)$. For the
last statement, $u_{0}$ satisfies the IVP

$u_{0}^{\prime}=-\frac{\sqrt{2}b}{\sqrt{5}}\cdot\sqrt{u_{0}^{5}+\gamma},\quad
u_{0}(0)=1,u_{0}(1)=0,$

with $\gamma>0$ is determined by the boundary conditions, and $\tilde{u}_{0+}$
satisfies

$u_{0+}^{\prime}=-\frac{\sqrt{2}b}{\sqrt{5}}u_{0+}^{5/2},\;u_{0+}(0)=1.$

Note that
$\frac{\sqrt{2}b}{\sqrt{5}}\cdot\sqrt{u^{5}+\gamma}-\frac{\sqrt{2}b}{\sqrt{5}}u^{5/2}$
is decreasing with $u$, hence the difference of slopes is the smallest at the
$x=0$ end and increases towards $x=1$ end of $[0,1]$. $\Box$

There is yet another way to numerically determine a better approximation for
$u_{0}$, by using the Runge-Kutta method and the estimates (3.3). For $t>0$ we
have $u_{\rm heat}(x)=t\cdot y_{\rm heat}(1-x)$, so we approximate ${\gamma}$
above by

  * $\displaystyle{\Gamma}:=\lim_{t\to 0^{+}}\biggl{(}t\cdot\frac{\sqrt{5}}{\sqrt{2}b}\tilde{y}_{+}^{\prime}(0)\biggr{)}^{2}={\biggl{(}1+\frac{1.5\sqrt{2}\,b}{\sqrt{5}}\biggr{)}}^{-10/3}=\biggl{[}\frac{\sqrt{5}}{\sqrt{2}\,b}\cdot(-\tilde{u}_{0+}^{\prime}(1))\biggr{]}^{2}.$ (6.1)

The numerical function $w_{0}$ defined by the IVP

  * $\displaystyle w_{0,\Gamma}^{\prime}=-\frac{\sqrt{2}\,b}{\sqrt{5}}\sqrt{w_{0,\Gamma}^{5}+\Gamma},\quad w_{0,\Gamma}(0)=1,$ (6.2)

satisfies the ODE $w_{0,\Gamma}^{\prime\prime}=b^{2}w_{0,\Gamma}^{4}$, the
same as $u_{0}$. It approximates the exact solution $u_{0}$ at least as well
as $\tilde{u}_{0+}$ does, and can be plotted by using the Runge-Kutta method.

Table 3. Numerical approximations of $u_{0}$ and their errors. | $b$ | $w_{0,\Gamma}(1)$ | $\tilde{u}_{0+}(1)-w_{0,\Gamma}(1)$
---|---|---
10 | 0.193 | 0.015
30 | 0.097 | 0.008
70 | 0.055 | 0.004
100 | 0.044 | 0.003
500 | 0.015 | 0.001
1000 | 9.5 E-3 | 7.8 E-4
50000 | 7.0 E-4 | 5.7 E-5
100000 | 4.4 E-4 | 3.6 E-5

numerical values for various $b$ | graphs of $w_{0,\Gamma},\;\tilde{u}_{0+}-w_{0,\Gamma},\;$ $b=10,100,1000,10000$

The table shows that
$\displaystyle\lim_{b\to\infty}(\tilde{u}_{0+}-w_{0,\Gamma})=0$. Nevertheless,
one may argue that, for low values of $b$, $w_{0,\Gamma}$ approximates loosely
the exact solution $u_{0}$. Fortunately, this numerical analysis also suggests
an easy way of obtaining truly sharp bounds upper and lower envelopes.

  * Procedure 6.2

Given a parameter $b$, one obtains sharp upper and lower envelopes by slightly
increasing the value of the exponent $-10/3=-3.33\dots$.

For an upper bound, we replace $-3.33$ by a slightly greater $r$ —let the
resulting function be $w_{0,r}$—, such that the value at $x=1$ is positive;
for lower bound, the value of $w_{0,r}$ at $x=0$ should be negative.

  * Example 6.3

(i) Let $b=10$. A few trials yield:

    * –

for $r_{-}=-2.136804$, we have $w_{0,r_{-}}(1)=-\,6.2\cdot 10^{-9}<0$;

    * –

for $r_{+}=-2.136805$, we have $w_{0,r_{+}}(1)=3.9\cdot 10^{-7}>0$.

Thus, for $b=10$, the exact solution $u_{0}$ is squeezed between the envelopes
$w_{0,r_{-}},w_{0,r_{+}}$, and the error is at most $4\cdot 10^{-7}$.

(ii) Let $b=10^{6}$.

    * –

$r_{-}=-3.129033$, $w_{0,r_{-}}=3.46\cdot 10^{-10}<0$;

    * –

$r_{+}=-3.129034$, $w_{0,r_{+}}=8.16\cdot 10^{-10}>0$.

For this choices, the exact solution is approximated with an error less than
$1.2\cdot 10^{-9}$.

(iii) Numerical data (see Table LABEL:tab:r) suggests that the correct value
of $r$ —the one separating the upper and the lower envelopes— follows the
rule:

    * $\displaystyle r\approx-\frac{10}{3}+\frac{\rho}{\ln(b)},\quad\rho\approx 2.8.$ (6.3)

So the exponent $-10/3$ in equation (6.1) defining $\Gamma$ is optimal, for
arbitrary $b$.

Table 4. Precision of the R-K method, for $r$ given by (6.1). $b$ | $\begin{array}[]{c}\rho_{-}=2.84\\\ w_{0,r_{-}}(1)\end{array}$ | $\begin{array}[]{c}\rho_{+}=2.8\\\ w_{0,r_{+}}(1)\end{array}$ | $\tilde{u}_{0+}(1)$
---|---|---|---
$10$ | -1.5 E-2 | -7.9 E-3 | 0.2087
$10^{2}$ | -8.6 E-5 | 1.4 E-3 | 0.0477
$10^{3}$ | -5.9 E-5 | 2.7 E-4 | 0.0103
$10^{4}$ | -2.2 E-5 | 4.9 E-5 | 2.2 E-3
$10^{5}$ | -6.0 E-6 | 9.5 E-6 | 4.8 E-4
$10^{6}$ | -1.4 E-6 | 1.8 E-6 | 1.0 E-4
$10^{7}$ | -3.4 E-7 | 3.7 E-7 | 2.2 E-5
$10^{8}$ | -7.9 E-8 | 7.6 E-8 | 4.8 E-6
$10^{9}$ | -1.7 E-8 | 1.5 E-8 | 1.0 E-6
$10^{10}$ | -3.9 E-9 | 3.2 E-9 | 2.2 E-7

The exact solution $u_{0}$ is squeezed between $w_{0,r_{-}}$ and $w_{0,r_{+}}$
(for $b\geqslant 16$). The value $b=10$ is an exception, one should take
$\rho_{+}=2.75$ rather than $2.8$. The table suggests that the error of the
computation is approximately $b^{-1}$.

##  7\. Numerical analysis

So far, we developed a theoretical framework. Here we confront our analysis
with numerical evidence, to probe our techniques.

###  7.1. Algorithm

The code below222It’s available also at:
https://drive.google.com/file/d/1sPbb11_43iZ6V_Q0vdWHw5h-LJkKo4Oz/view?usp=sharing.,
written in MAPLE, organizes the various formulae, and shows the
interdependence of the quantities introduced so far. It runs on a usual
computer, and allows the interested reader to experiment with the own
favourite parameters.

Note that if, for instance, one is interested in a partial lower bound for
$y\in[T^{3/4},T]$, it suffices to change $sss$ —first line of the code— to
$t^{3/4}$. The estimates will be more accurate. The value of $B$ is supposed
to be at least $1$.

Table 5. Code for upper/lower and global/partial envelopes
$\text{restart}:N:=20000:Digits:=N:\textbf{b := 70: t := 0.1: sss :=
$t^{1/2}$: }B:=bt^{1.5}/\sqrt{5}:T:=1/t:L:=1.5\sqrt{2}B:$
---
$BB:=bb*tt^{1.5}/\sqrt{5}:TT:=1/tt:tL:=L/qq:c0:=\text{write 30 iterations of
the $\tanh(LL*)$ function}:$
$ccm:=\tanh(LL*c0+\tanh^{-1}(c0/zz)):ccp:=tanh(LL*c0+\tanh^{-1}(1/zz)/(1-LL/\cosh(LL*c0)^{2})):$
$C:=2*cc/(1+ee-(1-ee)*cc):ZZ:=(T^{3/2}+((1-ee)*(1/2))*C)/(1+((1-ee)*(1/2))*C):$
$RR:=((1-(1/4)*q^{2}*(1+ee)^{2})*5)*y^{4}+7*CC*(1-ee)*y^{5/2}-2*CC^{2}*(-ee^{2}+4*ee-1)*y$
$-CC^{3}*ee*(1-ee)/y^{(}1/2)-CC^{4}*ee^{2}/y^{2}+q^{2}*(5*(1+ee)^{2}*(1/4)):$
$nom:=5*y^{13/2}*q^{2}+16*CC^{2}*y^{7/2}-5*q^{2}*y^{5/2}+14*CC*y^{5}+2*CC^{3}*y^{2}$
$-2*(70*CC*y^{23/2}*q^{2}+10*CC^{3}*y^{17/2}*q^{2}+36*CC^{3}*y^{17/2}-70*CC*y^{15/2}*q^{2}+36*CC^{5}*y^{11/2}$
$+25*y^{13}*q^{2}-10*CC^{3}*y^{9/2}*q^{2}+60*CC^{2}*y^{10}*q^{2}-5*CC^{4}*y^{7}*q^{2}+4*CC^{3}*y^{9/2}*r+9*CC^{2}*y^{10}$
$-5*y^{9}*q^{2}*r+54*CC^{4}*y^{7}-25*y^{9}*q^{2}+9*CC^{6}*y^{4}-60*CC^{2}*y^{6}*q^{2}+5*CC^{4}*y^{3}*q^{2}$
$+8*CC^{2}*y^{6}*r-4*CC^{4}*y^{3}*r+5*y^{5}*q^{2}*r)^{1/2}:$
$denom:=5*y^{13/2}*q^{2}-8*CC^{2}*y^{7/2}-5*q^{2}*y^{5/2}+4*CC^{4}*y^{1/2}-4*CC^{3}*y^{2}:$
$ey:=-nom/denom:$ # This is the root close to $1$ of the equation
$R(\dots)=r$:
$mxder:=(1-5*ss^{4}+(4+Dta)*ss^{5})^{1/2}:qqss:=2*(1+(1-ee)*CC*ss^{1.5}-ee*CC^{2}*ss^{3})/((1+ee)*mxder):$
$maxders:=subs(ss=sss,mxder):maxdert:=subs(ss=t,mxder):minder:=(1-5*t^{4}+4*t^{5})^{1/2}:$
$qmin:=1-(3/5)*t^{3}:emax:=1-(4/5)*t^{3}+((8-3*t)*(1/5))*t^{4}:emin:=1-(4/5)*t^{3}+(3/5)*t^{4}:$
#Define the global upper bound $typ$.
$tqp:=(((1-t^{3})*(t^{3}+5))/(5*(1-t^{4})))^{1/2}:tBp:=B/tqp:typ:=(\tanh(\tanh^{-1}(t^{1.5})+1.5*\sqrt{2}*tBp*(1-x)))^{-2/3}:$
#Define the global lower bound $tym$.
$emm:=0.73:cmin:=subs(LL=0.99*L,c0):C0mm:=1/emm:Dta:=\frac{((subs(x=0.25,typ)-1)}{(0.25*\sqrt{2}*B))^{2}}:$
$qqq:=subs(ee=emm,CC=C0mm,ss=t,qqss):qmm:=min(1.1,qqq):Lmm:=subs(qq=qmm,tL):$
$zmm:=subs(ee=emm,cc=cmin,ZZ):cmm:=evalf[N](subs(LL=Lmm,zz=zmm,ccm)):$
$Cmm:=evalf[N](subs(ee=emm,cc=cmm,C)):$
$\displaystyle
tym:=\Bigl{(}Cmm*\frac{T^{3/2}*(1+emm*\exp(2*Lmm*Cmm*(1-x)))+emm*Cmm*(\exp(2*Lmm*Cmm*(1-x))-1)}{T^{3/2}*(\exp(2*Lmm*Cmm*(1-x))-1)+Cmm*(\exp(2*Lmm*Cmm*(1-x))+emm)}\Bigr{)}^{2/3}:$
#Define the partial upper bound $yp$.
$\displaystyle
qqp:=\frac{2*(1+(1-ee)*CC*t^{1.5}-ee*CC^{2}*t^{3})}{(1+ee)*minder}:e1p:=subs(y=T,q=qmin,CC=cmin,r=0,ey):$
$q1p:=subs(ee=emin,CC=cmin,qqp):L1p:=subs(qq=q1p,tL):z1p:=subs(ee=e1p,cc=1,ZZ):$
$c1p:=evalf[N](subs(LL=L1p/e1p,zz=z1p,ccp)):C1p:=evalf[N](subs(ee=e1p,cc=c1p,C)):$
$R1:=evalf[N](subs(y=T,ee=e1p,q=q1p,CC=C1p,RR)):e2p:=e1p:q2p:=q1p+(1/15)*R1*t^{4}:$
$L2p:=subs(qq=q2p,tL):z2p:=subs(ee=e2p,cc=1,ZZ):c2p:=evalf[N](subs(LL=L2p/e2p,zz=z2p,ccp)):$
$C2p:=evalf[N](subs(ee=e2p,cc=c2p,C)):R2:=evalf[N](subs(y=T,ee=e2p,q=q2p,CC=C2p,RR)):$
$e3p:=e2p:q3p:=q2p+(1/15)*R2*t^{4}:L3p:=subs(qq=q3p,tL):z3p:=subs(ee=e3p,cc=1,ZZ):$
$c3p:=evalf[N](subs(LL=L3p/e3p,zz=z3p,ccp)):C3p:=evalf[N](subs(ee=e3p,cc=c3p,C)):$
$ep:=e3p:qp:=q3p:Cp:=C3p:Lp:=L3p:$
$\displaystyle
yp:=\Bigl{(}Cp*\frac{(T^{3/2}*(1+ep*\exp(2*Lp*Cp*(1-x)))+ep*Cp*(\exp(2*Lp*Cp*(1-x))-1))}{(T^{3/2}*(\exp(2*Lp*Cp*(1-x))-1)+Cp*(\exp(2*Lp*Cp*(1-x))+ep))}\Bigr{)}^{2/3};$
#Define the partial lower bound $yp$.
$qmax:=1.25:ccmin:=subs(LL=L/qmax,c0):CCmax:=1/emin:CCmin:=ccmin:qtm:=subs(ss=t,qqss):$
$qsm:=subs(ss=sss,qqss):ee0m:=emax:qq0m:=subs(ee=emax,CC=CCmax,qsm):$
$LL0m:=subs(qq=qq0m,tL):zz0m:=subs(ee=ee0m,cc=ccmin,ZZ):$
$cc0m:=evalf[N](subs(LL=LL0m,zz=zz0m,ccm)):CC0m:=evalf[N](subs(ee=ee0m,cc=cc0m,C)):$
$ee1m:=emax:qq1m:=subs(ee=emax,CC=CC0m,qsm):LL1m:=subs(qq=qq1m,tL):$
$zz1m:=subs(ee=ee1m,cc=ccmin,ZZ):cc1m:=evalf[N](subs(LL=LL1m,zz=zz1m,ccm)):$
$CC1m:=evalf[N](subs(ee=ee1m,cc=cc1m,C)):em:=ee1m:qqmm:=qq1m:Cm:=CC1m:Lm:=LL1m:$
$\displaystyle
ym:=\Bigl{(}Cm*\frac{T^{3/2}*(1+em*\exp(2*Lm*Cm*(1-x)))+em*Cm*(\exp(2*Lm*Cm*(1-x))-1)}{T^{3/2}*(\exp(2*Lm*Cm*(1-x))-1)+Cm*(\exp(2*Lm*Cm*(1-x))+em)}\Bigr{)}^{2/3}:$
$plot([typ-tym,yp-ym],x=0..1);plot([yp,ym],x=0..1);$

###  7.2. Numerical data

We are going to illustrate the enhanced precision of the partial bounds,
compared to the global ones. Below we plotted the graphs of various
differences between the global and the partial envelopes, and also the plots
of the envelopes themselves. The abbreviation B-L stands for ‘boundary layer’.

Table 6. Comparison of global and partial envelopes $\tilde{y}_{+}-\tilde{y}_{-}$ (red), $y_{+}-y_{-}$ (blue) |  $\tilde{y}_{+}-y_{+}$ (red), $y_{-}-\tilde{y}_{-}$ (blue) |  $\tilde{y}_{+}$ (red), $\tilde{y}_{-}$ (blue)
---|---|---
$\begin{array}[]{l}b=70\\\ t=0.1\\\ B=0.99\\\ bt^{0.25}=39\\\ \text{B-L no}\end{array}$ |  |
$\begin{array}[]{l}b=500\\\ t=0.1\\\ B=7.07\\\ bt^{0.25}=281\\\ \text{B-L yes}\end{array}$ |  |
$\begin{array}[]{l}b=10^{6}\\\ t=2\cdot 10^{-4}\\\ B=1.26\\\ bt^{0.25}=1.89\cdot 10^{5}\\\ \text{B-L yes}\end{array}$ |  |

It’s apparent that $\tilde{y}_{+}$ and $y_{+}$ are very close to each other
—in the second column of plots, the red graph basically are on the X-axis—,
while $y_{-}-\tilde{y}_{-}$ is not negligible; $y_{-}$ considerably improves
$\tilde{y}_{-}$. Moreover, the difference $y_{+}-y_{-}$ of the partial
envelopes —the blue graphs in the first column— is typically small. Thus the
work done before pays off. We included the example, $b=10^{6},t=2\cdot
10^{-4}$, to show that our methods is functional even in numerically
challenging situations.

The following table summarizes the computational error, the maximal difference
between the upper and lower envelopes, for various values of the parameters
$b,t$. The data suggests that the precision of the computation is related to
the value $bt^{0.25}$ of the boundary layer:

  * –

When computing with global envelopes, the error (for $y$) increases with it,
but the error for $u$ decreases.

  * –

When using partial envelopes, the precision increases (error decreases) with
the value of the boundary layer, both for $y$ and $u$. The numerical data
suggests that the number of exact decimals is of the same order as the
boundary-value $bt^{0.25}$. This is precisely the desired behaviour.

Table 7. Maximal computational error for $y_{\rm heat}$ and $u_{\rm heat}$. parameters and errors | max. difference | max difference
---|---|---
for $y$ and $u$ | global envelopes | partial envelopes
$b$ | $t$ | $B$ | $bt^{0.25}$ |  Xof $y$X | of $u$ | of $y$ | of $u$
500 | 0.1 | 7.07 | 281 | 0.16 | 0.016 | 0.027 | 0.0027
700 | 0.2 | 28 | 468 | 0.11 | 0.022 | 0.024 | 0.0048
5000 | 0.01 | 2.23 | 1581 | 1.1 | 0.011 | 0.017 | 1.7 E-4
10000 | 0.005 | 1.58 | 2659 | 2.17 | 0.01 | 0.013 | 6.5 E-5
1000000 | 0.0002 | 1.26 | 1.19 E+5 | 5.87 | 0.005 | 1.4 E-3 | 2.8 E-7

Here we listed the worst possible error. But, for a specific value $x$, the
computational error for $u_{\rm heat}(x)$ is $t\cdot(y_{+}(1-x)-y_{-}(1-x))$,
which can be much less than the maximal error, see Table 8. Moreover, one can
combine the envelope technique with the Runge-Kutta method in Section 3.4,
whose precision is $O(t\cdot\exp(-3B))$.

Table 8. Upper and lower bounds for $u_{\rm heat}$, for $b=500,t=0.1$. |  |  | number of
---|---|---|---
$x$ | $t\cdot y_{+}(x)$ | $t\cdot y_{-}(x)$ | exact decimals
$10^{-2}$ | 0.314022890404343 | 0.311867729652350 | 2
$10^{-4}$ | 0.969598013211494 | 0.969224267767387 | 3
$10^{-6}$ | 0.999684111921659 | 0.999680075693138 | 5
$10^{-8}$ | 0.999996839882428 | 0.999996799488355 | 6
$10^{-10}$ | 0.999999968398700 | 0.999999967994756 | 8

##  8\. Conclusion

The article investigates the equation of heat conduction with radiation:

$u^{\prime\prime}_{\rm heat}=b^{2}(u^{4}_{\rm heat}-t^{4}),\quad u_{\rm
heat}(0)=1,\;u_{\rm heat}(1)=t\in(0,1).$

It combines analytical with numerical techniques.

On the analytic side, we determined global and partial, upper and lower bounds
for the exact solution $u_{\rm heat}$. Also, we clarified the law, the value,
governing the boundary layer property of $u_{\rm heat}$. Moreover, we
estimated the initial derivative $u^{\prime}_{\rm heat}(1)$ of the solution to
the heat equation, which is essential for running the Runge-Kutta algorithm.
The precision of the estimate for $u_{\rm heat}$ is $O(t\cdot e^{-3B}).$

On the numerical side, we implemented the analytical arguments into the MAPLE
computer program and produced numerical data for various parameters. The
precision in estimating $u_{\rm heat}$ —the number of exact decimals— has the
same order of magnitude as the boundary-value $bt^{0.25}$.

## References

  * [1] Arpaci V., Conduction Heat Transfer. Addison-Wesley, Reading, Massachusetts, (1966).
  * [2] Carslaw H., Jaeger J., Conduction of Heat in Solids 2ed.. Oxford Science Publications, Clarendon Press, Oxford University Press, New York, (1988).
  * [3] Cannon J.R. The One-Dimensional Heat Equation. Encyclopedia of Math. and Appl. vol.23, Addison-Wesley, (1984).
  * [4] Firnett P., Troesch B, Shooting-splitting method for sensitive two-point boundary value problems. In: Proc. Conf. Numerical Solution of Ordinary Differential Equations, 19-20 October 1972, Lect. Notes Math. 362, Springer-Verlag, Berlin-Heidelberg-New York, 1974, pp. 408-433.
  * [5] Halic M., Tajarod R., Temimi H., An Accurate Analytical Solution to Troesch’s Problem Through Lower and Upper Envelope Techniques. Chaos, Solitons & Fractals: X 11 (2023), 100096.
  * [6] Holman J.P., Heat Transfer 10ed.. McGraw-Hill Series in Mech. Eng., McGraw-Hill, (2009).
  * [7] Mann W.R., Wolf F., Heat transfer between solids and gases under nonlinear boundary conditions. Quart. Appl. Math. 9 (1951), 163-184.
  * [8] Huang J., Cao L., Global regularity and multiscale approach for thermal radiation heat transfer. Multiscale Model. Simul. 12 (2014), 694-724.
  * [9] Modest M., Radiative Heat Transfer 3ed.. Academic Press, Elsevier New York London, (2013).
  * [10] Qin Y., Liu X., Yang X., Global existence and exponential stability for a 1D compressible and radiative MHD flow. J. Diff. Equations 253 (2012), 1439-1488.
  * [11] Tao L.N. Heat conduction with nonlinear boundary condition. J. Appl. Math. Phys. (ZAMP) 32 (1981), 144-155.
  * [12] Tokuda N., A new application of Lagrange-Bürmann expansions. II. Application to unsteady heat conduction problems with radiation. J. Appl. Math. Phys. (ZAMP) 34 (1983), 787-806.