Approximation of the modified error function

TL;DR

This paper develops explicit approximations for the modified error function related to heat transfer problems, demonstrating convergence to the classical error function as a key parameter approaches zero.

Contribution

It introduces a novel approximation method for the modified error function using power series in the parameter, applicable to phase-change heat transfer models.

Findings

Accurate approximations involving error and exponential functions.

Modified error function converges to the classical error function as delta approaches zero.

Characterizes the properties of the modified error function for small positive delta.

Abstract

In this article, we obtain explicit approximations of the modified error function introduced in Cho, Sunderland. Journal of Heat Transfer 96-2 (1974), 214-217, as part of a Stefan problem with a temperature-dependent thermal conductivity. This function depends on a parameter $\delta$, which is related to the thermal conductivity in the original phase-change process. We propose a method to obtain approximations, which is based on the assumption that the modified error function admits a power series representation in $\delta$. Accurate approximations are obtained through functions involving error and exponential functions only. For the special case in which $\delta$ assumes small positive values, we show that the modified error function presents some characteristic features of the classical error function, such as monotony, concavity, and boundedness. Moreover, we prove that the modified error function converges to the classical one when $\delta$ goes to zero.

Equations

$$\mathrm{erf}(x)=\frac{2}{\sqrt{\pi }}{\int }_0^x\exp (-{\xi }^2)d\xi x>0.$$

$$k(T)=k_0\left\{1+\delta \left(\frac{T-T_{\infty }}{T_m-T_{\infty }}\right)\right\},$$

$$[(1+\delta y(x))y^{\prime }(x){]}^{\prime }+2x{y}^{\prime }(x)=00<x<+\infty$$

$$y(0)=0$$

$$y(+\infty )=1.$$

$$\frac{x}{2}(1+x{)}^{3/2}(3+x)[1+(1+x{)}^{3/2}]=1x>0.$$

$${\Phi }_{\delta }(x)=\sum _{n=0}^{\infty }{\phi }_n(x){\delta }^{n}x>0,$$

$${\Psi }_{\delta ,m}(x)=\sum _{n=0}^m{\phi }_n(x){\delta }^{n}x>0\text{for }m\in {\mathbb{N}}_0.$$

$$\begin{array}{rl}& \sum _{n=0}^{\infty }\sum _{m=0}^{\infty }{\delta }^{n+m+1}\left({\phi }_n^{\prime }(x){\phi }_m^{\prime }(x)+{\phi }_n(x){\phi }_m^{\prime \prime }(x)\right)\\ & +\sum _{n=0}^{\infty }{\delta }^n\left({\phi }_n^{\prime \prime }(x)+2x{\phi }_n^{\prime }(x)\right)=0x>0.\end{array}$$

$$a(x,n,m)={\phi }_n^{\prime }(x){\phi }_m^{\prime }(x)+{\phi }_n(x){\phi }_m^{\prime \prime }(x)$$

$$b(x,n)={\phi }_n^{\prime \prime }(x)+2x{\phi }_n^{\prime }(x)$$

$$\sum _{n=1}^{\infty }\left(\sum _{k=1}^{n}a(x,k-1,n-k)+b(x,n)\right){\delta }^n+b(x,0)=0,x>0.$$

$$\sum _{k=1}^{n}a(x,k-1,n-k)+b(x,n)=0x>0,\forall n\in \mathbb{N},b(x,0)=0,x>0,$$

$${\phi }_n(0^{+})=0\forall n\in {\mathbb{N}}_0,{\phi }_0(+\infty )=1,{\phi }_n(+\infty )\forall n\in \mathbb{N}.$$

$$2x{\phi }_0^{\prime }(x)+{\phi }_0^{\prime \prime }(x)=0$$

$${\phi }_0(0^{+})=0$$

$${\phi }_0(+\infty )=1$$

$$2x{\phi }_n^{\prime }(x)+{\phi }_n^{\prime \prime }(x)=-\sum _{k=1}^n\left({\phi }_{k-1}^{\prime }(x){\phi }_{n-k}^{\prime }(x)+{\phi }_{k-1}(x){\phi }_{n-k}^{\prime \prime }(x)\right)$$

$${\phi }_n(0^{+})=0$$

$${\phi }_n(+\infty )=0,$$

$${\Psi }_{\delta ,0}(x)=\mathrm{erf}(x)x>0.$$

$${ϵ}_{\delta ,0}\to 0\text{and}\delta \to 0^{+},$$

$${ϵ}_{\delta ,0}=||{\Phi }_{\delta }-\mathrm{erf}|{|}_{\infty }.$$

$${\tau }_{\delta }(h)(x)=C_{\delta ,h}{\int }_0^x\frac{1}{1+\delta h(\eta )}\exp \left(-2{\int }_0^{\eta }\frac{\xi }{1+\delta h(\xi )}d\xi \right)d\eta x\ge 0,$$

$$C_{\delta ,h}={\left({\int }_0^{+\infty }\frac{1}{1+\delta h(\eta )}\exp \left(-2{\int }_0^{\eta }\frac{\xi }{1+\delta h(\xi )}d\xi \right)d\eta \right)}^{-1}.$$

$$\exists L>0/\forall {\delta }_1,{\delta }_2\in [0,{\delta }_0):||{\Phi }_{{\delta }_1}-{\Phi }_{{\delta }_2}|{|}_{\infty }\le L|{\delta }_1-{\delta }_2|,$$

$$x_1={\int }_0^{\eta }\frac{\xi }{1+{\delta }_{1}h(\xi )}d\xi ,x_2={\int }_0^{\eta }\frac{\xi }{1+{\delta }_{2}h(\xi )}d\xi (\eta >0\text{ fixed}).$$

$$\left|f(x_1)-f(x_2)\right|=|f^{\prime }(u)||x_1-x_2|,$$

$$|f^{\prime }(u)|\le |f^{\prime }(x_1)|\le 2\exp \left(-\frac{{\eta }^2}{1+{\delta }_0}\right)\text{ since }||h|{|}_{\infty }\le 1,$$

$$|x_1-x_2|\le \frac{{\eta }^2}{2}|{\delta }_1-{\delta }_2|.$$

$$\left|f(x_1)-f(x_2)\right|\le |{\delta }_1-{\delta }_2|{\eta }^2\exp \left(-\frac{{\eta }^2}{1+{\delta }_0}\right).$$

$$\begin{array}{rl}\left|\frac{f(x_1)}{1+{\delta }_{1}h(\eta )}-\frac{f(x_2)}{1+{\delta }_{2}h(\eta )}\right|& =\left|\frac{f(x_1)-f(x_2)}{1+{\delta }_{1}h(\eta )}+\frac{f(x_2)h(\eta )({\delta }_2-{\delta }_1)}{(1+{\delta }_{1}h(\eta ))(1+{\delta }_{2}h(\eta ))}\right|\\ & \le |f(x_1)-f(x_2)|+|f(x_2)||{\delta }_1-{\delta }_2|\\ & \le |{\delta }_1-{\delta }_2|\exp \left(-\frac{{\eta }^2}{1+{\delta }_0}\right)({\eta }^2+1).\end{array}$$

$$\begin{array}{rl}& \left|\frac{1}{C_{{\delta }_1,h_1}}-\frac{1}{C_{{\delta }_2,h_2}}\right|\le \left|\frac{1}{C_{{\delta }_1,h_1}}-\frac{1}{C_{{\delta }_1,h_2}}\right|\\ & \\ & +\left|{\text{bigintss}}_0^{+\infty }\frac{\exp \left(-2{\int }_0^{\eta }\frac{\xi }{1+{\delta }_1{h}_2(\xi )}d\xi \right)}{1+{\delta }_1{h}_2(\eta )}-\frac{\exp \left(-2{\int }_0^{\eta }\frac{\xi }{1+{\delta }_2{h}_2(\xi )}d\xi \right)}{1+{\delta }_2{h}_2(\eta )}d\eta \right|.\end{array}$$

$$\begin{array}{rl}|{\Phi }_{{\delta }_1}(x)-{\Phi }_{{\delta }_2}(x)|& \le \left|C_{{\delta }_1,{\Phi }_{{\delta }_1}}-C_{{\delta }_2,{\Phi }_{{\delta }_2}}\right|{\int }_0^{\infty }H({\Phi }_{{\delta }_2},{\delta }_2)(\eta )d\eta \\ & +C_{{\delta }_1,{\Phi }_{{\delta }_1}}\left|{\int }_0^{x}H({\Phi }_{{\delta }_1},{\delta }_1)(\eta )-H({\Phi }_{{\delta }_2},{\delta }_1)(\eta )d\eta \right|\\ & +C_{{\delta }_1,{\Phi }_{{\delta }_1}}\left|{\int }_0^{x}H({\Phi }_{{\delta }_2},{\delta }_1)(\eta )-H({\Phi }_{{\delta }_2},{\delta }_2)(\eta )d\eta \right|\forall x>0,\end{array}$$

$$H(h,\delta )(x)=\frac{\exp \left(-2{\int }_0^x\frac{\xi }{1+\delta h(\xi )}d\xi \right)}{1+\delta h(x)}x>0,$$

Full text

Keywords: Modified error function, error function, phase-change problem, temperature-

dependent thermal conductivity, nonlinear second order ordinary differential equation.

2010 AMS Subject Classification: 35R35, 80A22, 34B15, 34B08.

1 Introduction

Phase-change processes are present in a broad variety of natural, technological and industrial situations [1, 5, 9, 17, 24, 25, 28]. Modelling them properly is then crucial for understanding or predicting the evolution of many physical processes. One common assumption when modelling phase-change processes is to consider constant thermophysical properties. Nevertheless, it is known that certain materials present properties which seem to obey other laws. Recently, some models including variable latent heat, density, melting temperature or thermal conductivity have been proposed in [2, 4, 21, 22, 34, 36, 44].

In this sense, in 1974, Cho and Sunderland presented a similarity solution for a Stefan problem in which the thermal conductivity is a linear function of the temperature distribution [13]. It is well known that similarity solutions to Stefan problems with constant coefficients can be expressed in terms of the error function $\operatorname{erf}$,

[TABLE]

In contrast to this, the solution obtained by Cho and Sunderland involves another function, which they have called modified error function. It was defined as the solution to a nonlinear boundary value problem, and its existence and uniqueness was recently proved in [10] for thermal conductivities with moderate variations. In spite of the latter, the modified error function was widely used for solving diffusion problems [8, 14, 23, 29, 32, 37, 40], even before it was formally introduced by Cho and Sunderland in 1974 [15, 45].

When phase-change processes come from technological or industrial problems, not only appropriate models are required but also their solutions (or, at least, some properties of them). Sometimes, when explicit solutions are not known, models are solved through numerical methods which are tested with experimental data. When the latter are not available, one common practice is to test numerical methods by applying them to another problem whose explicit solution is known. Thus, having explicit solutions to models for phase-change processes is sometimes quite useful. Many works have been done in this direction, see for example [3, 6, 7, 11, 12, 16, 18, 19, 20, 26, 27, 30, 31, 35, 38, 39, 41, 42, 43, 44, 46]. Regarding the model in [13], explicit solutions are not known yet. Aiming to make a contribution in this sense, the main goal of this article is to propose some approximations of the modified error function.

In order to present our ideas clearly, we briefly recall how the modified error function arises from the original phase-change process. For simplicity, we consider the case of a one-phase melting problem for a semi-infinite slab with phase-change temperature $T_{m}$, whose boundary $x=0$ is maintained at a constant temperature $T_{\infty}>T_{m}$. For this case, the thermal conductivity from Cho and Sunderland is

[TABLE]

where $k_{0}>0$ is the thermal conductivity at $x=0$, and $\delta$ is some dimensionless parameter. Since $T=T_{m}$ at the free boundary, $\delta>-1$ becomes a necessary condition to assure the thermal conductivity is positive when $x=s(t)$. When the temperature distribution is assumed to be in the form $T(x,t)=A+B\Phi_{\delta}\left(\frac{x}{2\sqrt{\alpha_{0}t}}\right)$,111$\alpha_{0}$ is the coefficient of diffusion at $t=0$. for $A$ and $B$ constant, one obtains that $\Phi_{\delta}$ can be found by solving the following nonlinear boundary value problem (see details in [13]):

[TABLE]

The solution $\Phi_{\delta}$ to this problem is the already mentioned modified error function. Some plots for $\Phi_{\delta}$ are shown in Figure 1. They were obtained by numerically solving problem (3) for $\delta=-0.9,-0.5,0,0.5,1,2$. For the special case in which $\delta=0$, which corresponds to a constant thermal conductivity, one finds that the modified error function coincides with the classical one. The coincidence is stronger than that shown from the numerical computations, since it can be easily proved that the error function is the only solution to problem (3) when $\delta=0$. As $\delta$ moves away from zero, the modified error function differs more and more from the classical one. Nevertheless, both functions seem to share some properties (such as non-negativity, boundedness and rapid convergence to 1 when $x\to+\infty$). Moreover, when $\delta>0$, the modified error function seems to be increasing and concave, as the error function is. Observe that $-1<\delta<0$ is related to thermal conductivities that decrease when temperature increases (e.g. lead, methanol), whereas $\delta>0$ corresponds to thermal conductivities that increase as the temperature does (e.g. glycerin, mercury).

Finally, we recall that the existence and uniqueness of $\Phi_{\delta}$ in the set of non-negative bounded analytic functions was recently proved in [10] for small positive values of $\delta$ (i.e. for increasing thermal conductivities that present moderate variations with respect to their initial value). Moreover, an upper bound $\delta_{0}$ for the parameter $\delta$ was characterized as the unique positive solution to the equation:

[TABLE]

The approximations for the modified error function proposed in this article are based on the assumption that $\Phi_{\delta}$ admits a power series representation in the parameter $\delta$. More precisely, we assume that there exist functions $\varphi_{n}$ defined on $\mathbb{R}^{+}$ such that:

[TABLE]

and look for approximations $\Psi_{\delta,m}$ of the form

[TABLE]

The organization of the article is as follows. First (Section 2), we formally characterize each function $\varphi_{n}$ as the solution to a linear boundary value problem for a second order differential equation. The latter is homogeneous when $n=0$, but presents a non-zero source term dependent on $\varphi_{k}$ for $k=0,\dots,n-1$, when $n\in\mathbb{N}$. Then (Section 3), we present the zero order approximation $\Psi_{\delta,0}$. We find that it is given by the classical error function, and we prove that $\Phi_{\delta}$ uniformly converges to $\Psi_{\delta,0}=\operatorname{erf}$ when $\delta$ goes to zero. Since theoretical results are presented, the proof will be given only for those values of $\delta$ for which existence and uniqueness of the modified error function is known (i.e. we only consider $\delta\to 0^{+}$). After that (Section 4), we present the first and second order approximations $\Psi_{\delta,1}$ and $\Psi_{\delta,2}$. We obtain that $\Psi_{\delta,1}$ can be explicitly written in terms of the error and exponential functions only. By contrast, $\Psi_{\delta,2}$ involves some integrals whose values (explicit dependence on $x$) are not yet available in the literature. We analyze numerical errors between the approximations $\Psi_{\delta,1}$, $\Psi_{\delta,2}$, and the modified error function $\Phi_{\delta}$ when $\delta$ assumes small positive values, thus the existence and uniqueness of the modified error is assured. From them, we conclude that the first order approximation is better than the approximation of order two. Moreover, we find that it is also better than the approximation of zero order. We present some plots comparing $\Psi_{\delta,1}$ and $\Phi_{\delta}$ for some values of $\delta>-1$, which show good agreement. Finally (Section 5) we restrict the arguments again to small positive values of $\delta$ and prove that the modified and classical error functions share the properties of being increasing, concave and bounded functions.

2 Formal series representation of the modified error function

This Section is devoted to obtain a formal characterization of the coefficients $\varphi_{n}$ in the power series representation of the modified error function $\Phi_{\delta}$ given by (5).

Let $\delta>-1$ and $x>0$ be given. When $\Phi_{\delta}$ is defined by (5), formal computations from equation (3a) yield

[TABLE]

By introducing the following notation:

[TABLE]

equation (7) can be written as

[TABLE]

Therefore, the function $\Phi_{\delta}$ defined by (5) is a formal solution to problem (3) if and only if the functions $\varphi_{n}$, $n\in\mathbb{N}_{0}$, are such that

[TABLE]

That is, if and only if $\varphi_{0}$ and $\varphi_{n}$, $n\in\mathbb{N}$, are solutions to

[TABLE]

and

[TABLE]

respectively.

Remark 1**.**

We find from (12a) that functions $\varphi_{k}$ must be known for $k=0,\dots,n-1$ in order to find $\varphi_{n}$.

The next two Sections are dedicated to present and analyse the approximations $\Psi_{\delta,m}$ when $m=0,1,2$ and their coefficients are solutions to problems (11), (12) (see (6)). Each function $\Psi_{\delta,m}$ will be referred to as approximation of order $m$.

3 Approximation of order zero

The approximation of order zero is $\Psi_{\delta,0}=\varphi_{0}$, where $\varphi_{0}$ is a solution of problem (11). We note that the latter coincides with (3) when $\delta=0$. Thus, as it was already mentioned, its unique solution is the error function. Hence

[TABLE]

The remaining part of this Section is devoted to prove that the modified error function uniformly converges to the classical one, when the parameter $\delta$ goes to zero. We will restrict our analysis for those values of $\delta$ for which existence and uniqueness of the modified error function is known, i.e. to small positive values of $\delta$ [10]. Thus, our main goal will be to prove that

[TABLE]

where $\epsilon_{\delta,0}$ is the error between the classical and modified error functions, defined by

[TABLE]

The fact that the error function satisfies problem (3) when $\delta=0$, suggests to analyse the dependence of problem (3) on the parameter $\delta$. We begin by recalling the main result in [10]:

Theorem 3.1**.**

Let $\delta_{0}$ be the only solution to equation (4), and let $0\leq\delta<\delta_{0}$ be given. Then there exist a unique solution $\Phi_{\delta}$ to problem (3) in the set $K$ of all non-negative analytic functions in $\mathbb{R}_{0}^{+}$ which are bounded by 1. Moreover, $\Phi_{\delta}$ is given as the unique fixed point of the operator $\tau_{\delta}$ from $K$ to $K$ defined by

[TABLE]

for $h\in K$, with $C_{\delta,h}$ given by

[TABLE]

Proof.

See [10]. ∎

Remark 2**.**

In **[10]**, the case $\delta=0$ was not considered. This is because it corresponds to the classical case, in which the thermal conductivity is constant and the modified error function is the classical one. Nevertheless, besides this physical consideration, all theorems in **[10]** are still valid when $\delta=0$. 3. 2.

From (2), we see that $dk(T)/dT=\delta k_{0}/(T_{m}-T_{\infty})$. Thus, the bound on $\delta$ established by Theorem 3.1 to prove the existence and uniqueness of the modified error function $\Phi_{\delta}$ determines a necessary condition on the data of the associated Stefan problem to obtain similarity solutions. This condition establishes that the velocity of change of the thermal conductivity with respect to changes on the temperature distribution must be controlled by some multiple of the ratio between the thermal conductivity $k_{0}$ at $x=0$, and the difference between the phase-change and boundary temperatures, $T_{m}$ and $T_{\infty}$. In other words, that $s<\delta_{0}k_{0}/(T_{m}-T_{\infty})$, where $s$ is the slope in the linear dependence of $k$ on $T$.

Definition 3.1**.**

We say that problem (3) is Lipschitz continuous on the parameter $\delta$ if

[TABLE]

where $\Phi_{\delta_{1}}$, $\Phi_{\delta_{2}}$ are the only solutions in $K$ to problem (3) with parameters $\delta_{1}$, $\delta_{2}$, respectively.

Thus, if problem (3) is Lipschitz continuous on $\delta$, we find that the modified error function $\Phi_{\delta}$ converges uniformly on $x>0$ to the classical error function $\operatorname{erf}$, when $\delta\to 0^{+}$. In other words, that $\epsilon_{\delta,0}\to 0$ when $\delta\to 0^{+}$. Before proving the Lipschitz dependence of problem (3) on $\delta$, we introduce some preparatory results in the following:

Lemma 3.1**.**

Let $\delta_{1},\delta_{2}\in[0,\delta_{0})$, $h,h_{1},h_{2}\in K$ and $0\leq x\leq+\infty$ be given. The following estimations hold:

2. a)

$\displaystyle\bigintss_{0}^{x}\left|\frac{\exp\left(-2\displaystyle\int_{0}^{\eta}\frac{\xi}{1+\delta_{1}h(\xi)}d\xi\right)}{1+\delta_{1}h(\eta)}-\frac{\exp\left(-2\displaystyle\int_{0}^{\eta}\frac{\xi}{1+\delta_{2}h(\xi)}d\xi\right)}{1+\delta_{2}h(\eta)}\right|d\eta\leq$* $\displaystyle\dfrac{\sqrt{\pi}}{4}(1+\delta_{0})^{1/2}(3+\delta_{0})|\delta_{1}-\delta_{2}|$* 3. 4. b)

$\displaystyle\left|\frac{1}{C_{h_{1},\delta_{1}}}-\frac{1}{C_{h_{2},\delta_{2}}}\right|\leq\frac{\sqrt{\pi}}{4}(1+\delta_{0})^{1/2}(3+\delta_{0})\left(\delta_{0}||h_{1}-h_{2}||_{\infty}+|\delta_{1}-\delta_{2}|\right).$**

Proof.

2. a)

Let $f$ be the real function defined on $\mathbb{R}^{+}_{0}$ by $f(x)=\exp(-2x)$, and

[TABLE]

It follows from the Mean Value Theorem applied to function $f$ that

[TABLE]

where $u$ is a real number between $x_{1}$ and $x_{2}$. Without any lost of generality, we assume that $\delta_{1}\geq\delta_{2}$. Then, $x_{1}\leq x_{2}$ and we find

[TABLE]

Therefore,

[TABLE]

Then, we find

[TABLE]

The final bound is now obtained by integrating the last expression. 3. b)

First, we find the estimation

[TABLE]

Taking into consideration that the first term in the right hand side of (19) is bounded by $\frac{\sqrt{\pi}}{4}\delta_{0}(1+\delta_{0})^{1/2}(3+\delta_{0})||h_{1}-h_{2}||_{\infty}$ (see [10, Lemma 2.1]), the desired bound follows from (19) and item a).

∎

Theorem 3.2**.**

Problem (3) is Lipschitz continuous on the parameter $\delta$.

Proof.

Let $\delta_{1},\delta_{2}\in[0,\delta_{0})$ be given, and let $\Phi_{\delta_{1}},\Phi_{\delta_{2}}\in K$ be the solutions to problem (3) with parameters $\delta_{1}$, $\delta_{2}$, respectively.

Exploiting the fact that $\Phi_{\delta_{i}}$ is the fixed point of the operator $\tau_{\delta_{i}}$ defined by (16), $i=1,2$, we find

[TABLE]

where we have written

[TABLE]

with $h=\Phi_{\delta_{i}},\,\delta=\delta_{j},\,i,j=1,2.$

From the estimations

1. i)

$C_{\delta_{1},\Phi_{\delta_{1}}}\leq\frac{2\left(1+\delta_{0}\right)}{\sqrt{\pi}}$ 2. ii)

$\left|C_{\delta_{1},\Phi_{\delta_{1}}}-C_{\delta_{2},\Phi_{\delta_{2}}}\right|=\left|\frac{1}{C_{\delta_{1},\Phi_{\delta_{1}}}}-\frac{1}{C_{\delta_{2},\Phi_{\delta_{2}}}}\right|C_{\delta_{1},\Phi_{\delta_{1}}}C_{\delta_{2},\Phi_{\delta_{2}}}$ 3. $\leq\frac{1}{\sqrt{\pi}}(1+\delta_{0})^{5/2}(3+\delta_{0})\left(\delta_{0}||\Phi_{\delta_{1}}-\Phi_{\delta_{2}}||_{\infty}+|\delta_{1}-\delta_{2}|\right)$

(see Lemma 3.1) 4. iii)

$\left|\int_{0}^{x}H(\Phi_{\delta_{1}},\delta_{1})(\eta)-H(\Phi_{\delta_{2}},\delta_{1})(\eta)d\eta\right|\leq$ 5. $\displaystyle\frac{\sqrt{\pi}}{4}\delta_{0}(1+\delta_{0})^{1/2}(3+\delta_{0})||\Phi_{\delta_{1}}-\Phi_{\delta_{2}}||_{\infty}$

(see [10, Lemma 2.1]) 6. iv)

$\left|\int_{0}^{x}H(\Phi_{\delta_{2}},\delta_{1})(\eta)-H(\Phi_{\delta_{2}},\delta_{2})(\eta)d\eta\right|\leq\frac{\sqrt{\pi}}{4}(1+\delta_{0})^{1/2}(3+\delta_{0})|\delta_{1}-\delta_{2}|$

(see Lemma 3.1),

and (20), we obtain

[TABLE]

where $C=\delta_{0}(1+\delta_{0})^{3/2}(3+\delta_{0})$.

Since $\delta_{0}$ is the solution to equation (4), we find

[TABLE]

Thus, $0<C<1$. From this and (21), we obtain

[TABLE]

∎

In conclusion, we have found that the approximation of zero order $\Psi_{\delta,0}$ is the classical error function. Furthermore, for those values of $\delta$ for which existence and uniqueness of the modified error function $\Phi_{\delta}$ is known, we found that $\Phi_{\delta}$ uniformly converges to $\Psi_{\delta,0}=\operatorname{erf}$ on $x>0$. In particular, the latter suggest that $\Phi_{\delta}\simeq\Psi_{\delta,0}=\operatorname{erf}$ for small positive values of $\delta$, in agreement to Figure 1.

4 Approximations of order one and two

The approximations of order one and two are given by

[TABLE]

respectively, where $\varphi_{0}$ is the solution to problem (11) and $\varphi_{1}$, $\varphi_{2}$ are the solutions to problem (12) for $n=1$ and $n=2$. From Section 3, we know that $\varphi_{0}=\operatorname{erf}$. It enables us to define the source term in the differential equation of problem (12) for $n=1$. This problem can be explicitly solved and, from its solution, it can be also defined and solved problem (12) for $n=2$. We summarize these results in the following

Theorem 4.1**.**

2. a)

The only solution $\varphi_{1}$ to problem (12) for $n=1$ is given by

[TABLE] 3. b)

The only solution $\varphi_{2}$ to problem (12) for $n=2$ is given by

[TABLE]

where $g_{2}$, $\operatorname{erfc}$, $\operatorname{erfi}$ are the real functions defined in $\mathbb{R}^{+}$ by

[TABLE]

Proof.

It follows from standard results in ordinary differential equations (see, e.g. [33]). ∎

Plots for $\varphi_{n}$, $n=0,1,2$, are shown in Figure 2.

We will now investigate the relation between the approximations $\Psi_{\delta,1}$, $\Psi_{\delta,2}$, and the modified error function $\Phi_{\delta}$. The analysis will be again limited to those values of $\delta$ for which it is known the existence and uniqueness of $\Phi_{\delta}$, i.e. to $0<\delta<\delta_{0}\simeq 0.2037$ [10]. In contrast to the analysis presented in Section 3, investigations here will be based on numerical computations.

Numerical values for $\Phi_{\delta}$ were obtained by solving problem (3) through the routine bvodeS implemented in Scilab. The problem was solved for the domain $[0,10]$, by considering a uniform mesh $\mathcal{P}$ with step size $10^{-2}$.

Let $\mathcal{E}_{\delta,m}$ be the discrete error between $\Phi_{\delta}$ and $\Psi_{\delta,m}$, defined by

[TABLE]

Figure 3a shows some plots of $\mathcal{E}_{\delta,m}$ for $m=0,1,2$ and $\delta\in[0,0.2]\subset[0,\delta_{0})$. On one hand, we find that $\Psi_{\delta,1}$ and $\Psi_{\delta,2}$ are better approximations of $\Phi_{\delta}$ than $\Psi_{\delta,0}$. On the other hand, we also find that $\Psi_{\delta,1}$ is better than $\Psi_{\delta,2}$. This, together with the fact that $\Psi_{\delta,1}$ admits an explicit representation in terms of error and exponential functions only (in contrast to $\Psi_{\delta,2}$, which involves some integrals that can not be explicitly computed), turns $\Psi_{\delta,1}$ the best approximation among those proposed in this article. Figure 3b shows the comparison between the $\Psi_{\delta,1}$ and the modified error function $\Phi_{\delta}$, for $\delta=0.2$. Though we are not able to find a complete explanation of why $\Psi_{\delta,1}$ approximates better $\Phi_{\delta}$ than $\Psi_{\delta,2}$, we suggest that the numerical implementation of the integrals in the definition of $\varphi_{2}$ might be introducing non-negligible perturbations.

Finally, in Figure 4 we present plots for the modified error function $\Phi_{\delta}$ and the approximation of first order $\Psi_{\delta,1}$ for different values of $\delta>-1$. Even when they do not belong to the interval $[0,\delta_{0})$ over which theoretical results on existence and uniqueness of $\Phi_{\delta}$ are known, very good agreement is obtained. This enforces the former conclusions of being $\Psi_{\delta,1}$ the best approximation of $\Phi_{\delta}$, among $\Psi_{\delta,m}$ for $m=0,1,2$.

5 Properties of the modified error function

We end this article by proving that the modified error function $\Phi_{\delta}$ found in [10] shares some basic properties with the classical error function $\operatorname{erf}$. More precisely, those of being an increasing concave non-negative and bounded function.

Theorem 5.1**.**

If $0<\delta<\delta_{0}$, then the only solution $\Phi_{\delta}$ in $K$ to problem (3) satisfies the following properties:

[TABLE]

Proof.

The first property in (26) is a direct consequence of the fact that $\Phi_{\delta}$ belongs to $K$. In order to prove the second one, we start by showing that $\Phi_{\delta}^{\prime}(x)\neq 0$ for all $x>0$. We will assume that there exists $x_{0}>0$ such that $\Phi_{\delta}^{\prime}(x_{0})=0$ and we will reach a contradiction. Since equation (3a) can be written as

[TABLE]

and we know that

[TABLE]

we find that $\Phi_{\delta}^{\prime\prime}(x_{0})=0$. From this, by differentiating (27) and taking (28) into consideration, it follows that $\Phi_{\delta}^{\prime\prime\prime}(x_{0})=0$. We continue in this fashion obtaining that $\Phi_{\delta}^{(n)}(x_{0})=0$ for all $n\in\mathbb{N}$. This implies that $\Phi_{\delta}\equiv 0$ in $\mathbb{R}_{0}^{+}$, since $\Phi_{\delta}$ is an analytic function. But the latter contradicts that $\Phi_{\delta}(+\infty)=1$. Therefore, $\Phi_{\delta}^{\prime}(x)\neq 0$ for all $x>0$. This implies that the function $\Phi_{\delta}^{\prime}$ does not change its sign in $\mathbb{R}_{0}^{+}$. Since $\Phi_{\delta}(0)\leq 1$ and $\Phi_{\delta}(+\infty)=1$, it follows that $\Phi_{\delta}^{\prime}(x)>0$ for all $x>0$. Finally, the last property in (26) follows straightforward from the previous ones and the fact that $\Phi_{\delta}^{\prime\prime}$ is given by

[TABLE]

∎

Remark 3**.**

Note that Figure 4a suggest that the assumption $\delta>0$ can not be removed from the statement of Theorem 5.1 since the two last properties in (26) seem to fail.

6 Conclusions

In this article, we have proposed a method to obtain approximations of the modified error function $\Phi_{\delta}$ introduced by Cho and Sunderland in 1974 [13] as part of a Stefan problem with variable thermal conductivity. This is defined as the solution to a nonlinear boundary value problem for a second order ordinary differential equation which depends on a parameter $\delta>-1$. By assuming that $\Phi_{\delta}$ admits a power series representation in $\delta$, we proposed some approximations $\Psi_{\delta,m}$ given as the partial sums of the first $m$ terms. It was presented three of them: the zeroth order approximation $\Psi_{\delta,0}=\operatorname{erf}$; the first order approximation $\Psi_{\delta,1}$, which can be written in terms of error and exponential functions only; and the second order approximation $\Psi_{\delta,2}$, which can be written in terms of the error and exponential functions, and some integrals of combinations of them. Analysis of errors between the approximations and the original function was performed by considering only those values of $\delta$ for which existence and uniqueness of the modified error function is known, i.e. to small positive values of the parameter $\delta$. When $m=0$ it was found that $\Phi_{\delta}$ uniformly converges to $\Psi_{\delta,0}$. This suggest that $\Phi_{\delta}\simeq\Psi_{\delta,0}$ for small values of $\delta$. When $m=1,2$, numerical investigations suggest that $\Psi_{\delta,1}$ and $\Psi_{\delta,2}$ are also accurate approximations for $\Phi_{\delta}$. In particular, it was obtained that $\Psi_{\delta,1}$ and $\Psi_{\delta,2}$ are better approximations than $\Psi_{\delta,0}$, and that $\Psi_{\delta,1}$ is better than $\Psi_{\delta,2}$. This, together with the simple expression of $\Psi_{\delta,1}$ in terms of the error and exponential functions only, turns the first order approximation the best one among those presented here. The fact that $\Psi_{\delta,1}$ seems to be a better approximation than $\Psi_{\delta,2}$ can not be completely addressed by the authors. Nevertheless, we suggest that the numerical implementation of the integrals in the definition of the second order approximation might be introducing non-negligible perturbations. Comparisons between $\Psi_{\delta,1}$ and $\Phi_{\delta}$ were also presented for values of $\delta>-1$. Good agreement was obtained, even for those values of $\delta$ which do not belong to the theoretical interval determined by the existence and uniqueness results of $\Phi_{\delta}$. Finally, we proved that the modified error function is an increasing non-negative concave function which is bounded by 1, as the classical error function is, provided $\delta$ assumes small positive values. Results presented here can be used to obtain explicit approximate solutions to Stefan problems for phase-change processes with linearly temperature-dependent thermal conductivity. This investigation suggest that the proof of existence and uniqueness of $\Phi_{\delta}$ for $-1<\delta<0$ and evaluation of integrals involving exp-, erf-, erfc-functions are still an open problem in the mathematical analysis of Stefan problems.

Acknowledgements

This paper has been partially sponsored by the Project PIP No. 0275 from CONICET-UA (Rosario, Argentina) and AFOSR-SOARD Grant FA 9550-14-1-0122. The authors would like to thanks the anonymous referees whose insightful comments have benefited the presentation of this article.