Solution of parabolic free boundary problems using transmuted heat polynomials

TL;DR

This paper introduces a numerical method for solving parabolic free boundary problems by utilizing transmutation operators to construct solution systems that generalize heat polynomials, with an implementation algorithm provided.

Contribution

The paper presents a novel numerical approach based on transmutation operators to efficiently solve free boundary problems for parabolic equations, extending heat polynomial methods.

Findings

Efficient construction of solution systems using transmutation operators.

Generalization of heat polynomials for free boundary problems.

Implementation algorithm demonstrated for the proposed method.

Abstract

A numerical method for free boundary problems for the equation \[ u_{xx}-q(x)u=u_t \] is proposed. The method is based on recent results from transmutation operators theory allowing one to construct efficiently a complete system of solutions for this equation generalizing the system of heat polynomials. The corresponding implementation algorithm is presented.

Equations

$$u_{xx}-q(x)u=u_t$$

$$u_{xx}-q(x)u=u_t.$$

$$\mathbf{A}=\frac{{\partial }^2}{\partial x^2}-q(x)$$

$$D(s)=\left\{(x,t)\in {\mathbb{R}}^2:0<x<s(t),0<t<T\right\},$$

$$\mathbf{A}u(x,t)=u_t(x,t),(x,t)\in D(s),$$

$$\left[{\mathbf{G}}_{1}u\right](x,0)$$

$$\left[{\mathbf{G}}_{2}u\right](0,t)$$

$$\displaystyle u\bigl{(}s(t),t\bigr{)}$$

$$\displaystyle u_{x}\bigl{(}s(t),t\bigr{)}$$

$$h_n(x,t)=\sum _{k=0}^{\left[n/2\right]}{c}_k^n{x}^{n-2k}{t}^k,\text{where }{c}_k^n=\frac{n!}{(n-2k)!k!}$$

$$h_0(x,t)$$

$$h_3(x,t)$$

$$u_{xx}=u_t$$

$$\left(\frac{d^2}{d{x}^2}-q(x)\right)f(x)=0,x\in (0,L),$$

$$f(0)=1.$$

$$X^{(0)}(x)\equiv 1,X^{(n)}(x)=n{\int }_0^x{X}^{(n-1)}(s){\left(f^2(s)\right)}^{(-1{)}^n}\mathrm{d}s,n=1,2,\dots$$

$${\tilde{X}}^{(0)}\equiv 1,{\tilde{X}}^{(n)}(x)=n{\int }_0^x{\tilde{X}}^{(n-1)}(s){\left(f^2(s)\right)}^{(-1{)}^{n-1}}\mathrm{d}s,n=1,2,\dots .$$

$${\phi }_n(x)=\{\begin{array}{ll}f(x)X^{(n)}(x),& n\text{ odd},\\ f(x){\tilde{X}}^{(n)}(x),& n\text{ even},\end{array}$$

$$\mathbf{T}u(x)=u(x)+{\int }_{-x}^{x}K(x,y)u(y)dy$$

$$\mathbf{AT}[u]=\mathbf{T}\left[\frac{d^{2}u}{d{x}^2}\right]$$

$$\mathbf{T}\left[1\right]=f.$$

$$\mathbf{T}\left[x^n\right]={\phi }_n(x)\text{for any }n\in {\mathbb{N}}_0.$$

$$\mathbf{A}u=u_t.$$

$$\underset{(x,t)\in D_c}{\max }\left|u(x,t)-\sum _{n=0}^N{a}_n{H}_n(x,t)\right|<\epsilon .$$

$$H_n(x,t)=\sum _{k=0}^{\left[n/2\right]}{c}_k^n{\phi }_{n-2k}(x)t^k.$$

$$u_N(x,t)=\sum _{n=0}^N{a}_n{H}_n(x,t)$$

$$s_K(t)=\sum _{k=0}^K{b}_k{\beta }_k(t).$$

$$v\left(\bar{y}\right)=\bigl{(}v\left(y_{0}\right),...,v\left(y_{N_{y}}\right)\bigr{)}^{T}.$$

$${\Vert v(\bar{y})\Vert }^2=\sum _{i=0}^{N_y}|v(y_i){|}^2.$$

$$I_1\left(\bar{a},\bar{b}\right)$$

$$I_2\left(\bar{a},\bar{b}\right)$$

$$I_3\left(\bar{a},\bar{b}\right)$$

$$I_4\left(\bar{a},\bar{b}\right)$$

$$F\left(\bar{a},\bar{b}\right)=\sum _{i=1}^4{I}_i^2\left(\bar{a},\bar{b}\right),$$

Full text

[TABLE]

is proposed. The method is based on recent results from transmutation operators theory allowing one to construct efficiently a complete system of solutions for equation (1) generalizing the system of heat polynomials. The corresponding implementation algorithm is presented.

1 Introduction

Free boundary problems (FBPs) for parabolic equations are of considerable interest in physics (e.g., the Stefan problem) and in financial mathematics (e.g., the problem of pricing of an American option). One of the relatively simple and practical methods proposed for solving FBPs involving the heat equation is the heat polynomials method (see [3], [31], [4], [30], [35], [11]) based on the fact that the system of the heat polynomials represents a complete family of solutions of the heat equation. In the book [3] D. Colton proposed to extend this method onto parabolic equations with variable coefficients by constructing an appropriate transmutation operator and obtaining with its aid the corresponding transmuted heat polynomials. However, the construction of the transmutation operator is a difficult problem itself. In the present work we show that the transmuted heat polynomials required for Colton’s approach can be constructed without knowledge of the transmutation operator by a simple and robust recursive integration procedure. To this aim a recent result from [2] concerning a mapping property of the transmutation operators is used.

This makes possible to extend the heat polynomials method onto equations of the form

[TABLE]

We mention that linear parabolic equations of a more general form with coefficients depending on one variable reduce to (2) (see, e.g., [3, Chap. 2] or [27]).

Thus, we propose a numerical method for approximate solution of a class of FBPs involving (2). This method of transmuted heat polynomials will be designated by THP. The main aim of this paper is to explain it in detail and to propose a simple to implement algorithm for its application.

The subject of FBPs appears in many different fields and applications, as such, presents a large variety of formulations and still open questions. We will not be focusing on the existence and uniqueness of the solution in this paper assuming that it exists in the classical sense.

The paper is structured as follows. In Section 2 we state the FBP. In Section 3 we define transmutation operators and present some of their properties motivating the development of the THP method. In Section 4 we construct the conceptual algorithm for the implementation of the THP method for solving FBPs. In Section 5 we consider an example with an exact solution, with its aid we illustrate the performance of the method. In Section 6 we discuss the possibilities of generalization of our construction and its application to more general FBPs.

2 Statement of the problem

Consider the differential expression

[TABLE]

with $q\in C[0,L]$, $q:[0,L]\rightarrow\mathbb{C}$.

Every $s\in C^{1}[0,T]$, such that $0<s(x)\leq L$ for all $x\in[0,T)$ and $s(0)=l$, defines a domain

[TABLE]

see Figure 1. Consider the first order linear differential operators $\mathbf{G}_{1}=\gamma_{11}(x)+\gamma_{12}(x)\frac{\partial}{\partial x}$ and $\mathbf{G}_{2}=\gamma_{21}(t)+\gamma_{22}(t)\frac{\partial}{\partial x}$ , where $\gamma_{ij}$ are some given continuous functions.

Problem 1

Find functions $(u,s)$ such that

(i)

$s\in C^{1}[0,T]$* and such that $0<s(x)\leq L$ for all $x\in[0,T)$ and $s(0)=l$,*

(ii)

u\in C^{2,1}\bigl{(}D(s)\bigr{)},

(iii)

the following equation is satisfied on $D(s)$

[TABLE]

(iv)

and the following boundary conditions are satisfied

[TABLE]

here the dot over the function $s$ means the derivative with respect to the variable $t$. The last condition is usually known as the equation of heat balance or as the Stefan condition.

This problem is broadly studied in literature (see, e.g., [8], [34], [6], [26], [7] and [36] for additional bibliography). In particular, the classical one dimensional one phase Stefan problem is a special case of Problem 1. Since the subject of this paper is the approximate numerical method for solution of Problem 1, we make the following assumption.

Assumption 2

There exists a unique solution to Problem 1.

A relevant example of an existence and uniqueness result is given in [7, Theorem 1], see Remark 11 below.

3 Transmuted heat polynomials

3.1 Heat polynomials

The heat polynomials are defined for $n\in\mathbb{N}_{0}$ as (see, e.g., [32] and [37])

[TABLE]

and $\left[\cdot\right]$ denotes the entire part of the number. The first five heat polynomials are

[TABLE]

The set of the heat polynomials $\left\{h_{n}\right\}_{n\in\mathbb{N}_{0}}$ represents a complete system of solutions for the heat equation

[TABLE]

on any domain $D(s)$ defined by (3) (see [5]).

3.2 Formal powers and the transmutation operator

Let $f$ be a nonvanishing (in general, complex valued) solution of the equation

[TABLE]

such that

[TABLE]

The existence of such solution111In fact the only reason for the requirement of the absence of zeros of the function $f$ is to make sure that the auxiliary functions (11) be well defined. As was shown in [20] this can be done even without such requirement, but corresponding formulas are somewhat more complicated. for any complex valued $q\in C[0,L]$ was proved in [19] (see also [1]).

Consider two sequences of recursive integrals (see [15], [16], [19])

[TABLE]

and

[TABLE]

Definition 3

The family of functions $\left\{\varphi_{n}\right\}_{n=0}^{\infty}$ constructed according to the rule

[TABLE]

is called the system of formal powers associated with $f$.

The formal powers arise in the spectral parameter power series (SPPS) representation for solutions of the Sturm-Liouville equation, see [14], [19], [12], [16].

The following result from [25] (see also [17] and [21] for additional details) and from [2] guarantees the existence of a transmutation operator associated with $f$ and shows its connection with the system of formal powers.

Theorem 4

Let $q\in C[0,L]$. Then there exists a unique complex valued function K(x,y)\in C^{1}\bigl{(}[0,L]\times[-L,L]\bigr{)} such that the Volterra integral operator

[TABLE]

defined on $C[0,L]$ satisfies the equality

[TABLE]

for any $u\in C^{2}[0,L]$ and

[TABLE]

Note that if $v=\mathbf{T}u$ then $v(0)=u(0)$ and $v^{\prime}(0)=u^{\prime}(0)+f^{\prime}(0)u(0)$.

Theorem 5** ([2])**

[TABLE]

3.3 Construction of the transmuted heat polynomials

Denote $H_{n}=\mathbf{T}\left[h_{n}\right]$. Thus, $H_{n}(x,t)=h_{n}(x,t)+\int_{-x}^{x}K(x,y)h_{n}(y,t)\,dy$. The functions $H_{n}$ are called the transmuted heat polynomials [18].

We have $\left(\mathbf{A}-\frac{\partial}{\partial t}\right)H_{n}=\mathbf{AT}h_{n}-\frac{\partial}{\partial t}\mathbf{T}h_{n}=\mathbf{T}\left(\frac{\partial^{2}}{\partial x^{2}}-\frac{\partial}{\partial t}\right)h_{n}=0$ and hence every $H_{n}$ is a solution of the equation

[TABLE]

The set $\left\{H_{n}\right\}_{n\in\mathbb{N}_{0}}$ is a complete system of solutions of (12) in the following sense.

Proposition 6

Let $u$ be a classical solution of (12) in $D(s)$, continuous in the closure $\bar{D}(s)$. Then for any compact set $D_{c}\subset D(s)$ and any given $\varepsilon>0$ there exist $N$ and constants $a_{0},\ldots,a_{N}$ such that

[TABLE]

Moreover, if $s(t)$ can be extended to a function analytic in the disk $\{t\in\mathbb{C}:|t|\leq T\}$, the uniform approximation property is valid on the whole set $\bar{D}$.

Proof. Suppose first that $q\in C^{1}[0,L]$. Note that any compact set $D_{c}$ can be covered by a subset $D_{a}\subset\bar{D}(s)$ of the form $D_{a}=\{(x,t)\in\mathbb{R}^{2}:0\leq x\leq x_{1}(t)\leq s(t),\ 0\leq t\leq T\}$ with analytic function $x_{1}(t)$. For the set $D_{a}$ the proof is completely similar to that of [3, Theorem 2.3.3] with the only change that the transmutation operator $\mathbf{T}$ and its inverse are used. For the general case $q\in C[0,L]$ one approximates $q$ by a $C^{1}$ function.

Theorem 7

The transmuted heat polynomials admit the following form

[TABLE]

Proof. This equality is an immediate corollary of Theorem 5. Indeed, we have $H_{n}(x,t)=\mathbf{T}h_{n}(x,t)=\sum_{k=0}^{\left[n/2\right]}c_{k}^{n}\mathbf{T}\left[x^{n-2k}\right]t^{k}$, where Theorem 5 is used.

The explicit form of the functions $H_{n}$ presented in this theorem makes possible the construction of the approximate solution to Problem 1 by the THP method.

4 Description of the method

We proceed to the step by step construction of the THP method, summarizing the algorithm at the end of the section.

Assume that we have already calculated the formal powers $\varphi_{k}$ and the functions $H_{n}$. Further, let $N$ be the highest index of the formal power considered, or equivalently the highest degree in $x$ of the considered heat polynomials. We denote by

[TABLE]

the approximation of the solution $u(x,t)$, and by $\bar{a}=\left(a_{0},...,a_{N}\right)^{T}$ the column-vector of the unknown coefficients. Note that by construction functions $u_{N}$ satisfy equation (4), all we need is to find suitable coefficients $\bar{a}$.

Denote by $\left\{t_{i}\right\}$ an ordered set of $N_{t}+1$ points of the interval $[0,T]$, with $t_{0}=0<t_{1}<...<t_{N_{t}}=T$, $\bar{t}=\left(t_{0},...,t_{N_{t}}\right)^{T}$, and by $\left\{x_{i}\right\}$ an ordered set of $N_{x}+1$ points of the initial boundary, with $x_{0}=0<x_{1}<...<x_{N_{x}}=l$, $\bar{x}=\left(x_{0},...,x_{N_{x}}\right)^{T}$.

Consider a set of $K+1$ linearly independent differentiable functions, $\beta_{k}:[0,T]\rightarrow\mathbb{R}$. We are looking for the free boundary in the form

[TABLE]

Denote the vector of the unknown coefficients222It is possible to search for the free boundary in a more general form, see Section 6.3. by $\bar{b}=\left(b_{0},...,b_{K}\right)^{T}$. For any function $v$ defined on the set of points $\left\{y_{i}\right\}_{i=0,1,...,N_{y}}$ the following notation is introduced

[TABLE]

A numerical approximation problem consists in finding a set of the coefficients $\left(\bar{a},\bar{b}\right)$ that best fits the conditions of Problem 1 with the following norm chosen

[TABLE]

The following magnitudes are to be minimized,

[TABLE]

Each of them is related to a boundary condition from (5)–(8). With introduction of the value function

[TABLE]

the minimization problem can be stated as follows.

Problem 8

*Find333For a function $f:X\rightarrow Y$, the $\arg\min$ over a subset $S$ of $X$ is defined as

$\underset{x\in S\subseteq X}{\arg\min}\,f(x):=\left\{x:x\in S\wedge\forall y\in S:f(y)\geq f(x)\right\}.$

[TABLE]

subject to

[TABLE]

Note that instead of the uniform norm chosen in [31] for the heat polynomials method, we used the $L_{2}$ norm in the value function (14). The main reason for such choice was to take advantage of the presence of the so-called separable linear parameters (see [33, Chap. 6.2], see also [10]) and to reduce the number of parameters in the value function. Indeed, for each fixed $\bar{b}$, the constrained minimization Problem 8 reduces to an unconstrained linear least squares problem for parameters $\bar{a}$ and can be easily solved exactly (see Subsection 4.1). That is, for each $\bar{b}$ we can define

[TABLE]

So instead of minimizing the value function $F$ over an $N+K+2$ dimensional space of parameters $(\bar{a},\bar{b})$, the problem can be reduced to minimization of the function

[TABLE]

over a $K+1$ dimensional space. Such reformulation of the original optimization problem leads to a more robust convergence of numerical optimization algorithms, c.f., [10]. We state this new minimization problem as follows.

Problem 9

Find

[TABLE]

subject to

[TABLE]

4.1 Linear minimization problem

Let a vector $\tilde{b}$ be fixed. We have to solve the linear least squares problem (16). Let us denote the corresponding approximate boundary by $\tilde{s}$,

[TABLE]

Note that problem (16) is equivalent to solving the overdetermined system of linear equations (resulting from the boundary conditions (5)–(8))

[TABLE]

where

[TABLE]

or in the matrix form

[TABLE]

where

[TABLE]

Note that the derivatives $\varphi^{\prime}$ in (20) do not require numerical differentiation and can be obtained in a closed form from (11).

The solution of this overdetermined system coincides [24, Thm. 5.14] with the unique solution of the following fully determined system of linear equations

[TABLE]

See also [23] or [28] for various methods of efficient solution of equation (21) or (22). In particular, we used the Moore-Penrose pseudo-inverse method.

4.2 Implementation

Here we present the algorithm of the implementation of the THP method for Problem 1.

Find a particular solution of equation (9) satisfying (10). For example, the SPPS method presented in [19] can be used. 2. 2.

Construct the formal powers on an interval from [math] to some $L>l$. In this paper we represented all the functions involved by their values on a uniform mesh and used the Newton-Cotes six point integration rule. Since we may need the values of the formal powers at arbitrary points $\tilde{s}\left(\bar{t}\right)\subset[0,L]$, we approximated the formal powers by splines passing through their values on the selected mesh. Spline integration can be used as well for the construction of the formal powers, c.f., [13]. See also [17] for the discussion of other possible methods. 3. 3.

Compute the coefficients $B_{n}$ and $C_{n}$ given by (18) and (19) and the respective functions $g_{1}$ and $g_{2}$. Since these conditions are independent of the free boundary (i.e. of the choice of $\tilde{b}$), they need to be computed only once. 4. 4.

Construct a function that solves problem (16) for each given set of coefficients $\tilde{b}$. 5. 5.

Construct the value function $\tilde{F}\left(b\right)$ using (14). 6. 6.

Solve the constrained minimization Problem 9 using any suitable algorithm, see, e.g., [28]. For the numerical illustration we used the Matlab function fmincon. Note that for faster convergence we may initially solve the minimization Problem 9 for some small value $K^{\prime}<K$ and use the coefficients $b^{0}_{k}$, $k=0,\ldots,K^{\prime}$ obtained as an initial value for the optimization algorithm for larger value $K$.

5 Numerical illustration

5.1 Example with an exact solution and a tractable free boundary

For the numerical experiment we constructed an example of an FBP admitting an exact solution.

For every $s\in C^{1}[0,1]$, with $s(0)=1$ consider the domain

[TABLE]

Problem 10

Find the pair $(s,u)$ such that $s\in C^{1}[0,1]$, $u\in C^{2,1}\left(D(s)\right)$ and

[TABLE]

with the following boundary conditions

[TABLE]

where

[TABLE]

and $\operatorname{Ei}^{-1}$ stands for the inverse function of the exponential integral $\operatorname{Ei}$,

[TABLE]

Problem 10 possesses an exact solution $u$ given by

[TABLE]

and the free boundary

[TABLE]

Remark 11

The existence and the uniqueness of the solution for $t\in\left[0,1\right]$ is guaranteed by [7, Theorem 1]. For the application of the theorem it is necessary to use the transformation $v(x,t)=u(x,t)+\gamma(t)$, where $\gamma(t)=\exp\left(-\operatorname{Ei}^{-1}\left(2C-2e^{-t}\right)-t\right)$.

5.2 Numerical illustration

We proceed by presenting numerical results delivered by the algorithm described in Section 4. All the calculations were carried out in Matlab R2012a.

A particular solution to equation (9) and the formal powers (11) were represented by their values on a 2000-points uniform mesh. The Newton-Cotes integration rule was used for their computation. The number of formal powers considered was set at $N=12$. Finally the function spapi was used to create splines approximating the formal powers. For computing the $\operatorname{Ei}^{-1}$ we have used a polynomial interpolation on a fine grid for the function $\operatorname{Ei}$. The values for the function $\operatorname{Ei}$ were constructed by the series expansions presented in [9, (8.214)], see also [29].

The sets of points $\left\{t_{i}\right\}$ and $\left\{x_{i}\right\}$ considered are both equidistant grids of $101$ points. The free boundary $s$ is sought in terms of polynomials. The condition $s(0)=1$ inspires the following form

[TABLE]

In our calculations we have considered $K=6$.

The linear problem (16) was solved by using the Matlab function pinv. We have found the solution to Problem 9 using the Matlab routine fmincon.

The proposed algorithm converged rapidly for various initial free boundaries tested. On Figure 2 we present one of the initial boundaries, the exact free boundary $s(t)$ calculated from equation (28) and the difference between the exact free boundary $s(t)$ and the obtained approximation $s_{K}(t)$. The obtained approximate boundary was

[TABLE]

here and after the coefficients are presented up to the eighth decimal place. The coefficients of the vector $\bar{a}$ that corresponds to the solution of the linear problem (16) with $\bar{b}$ consisting of the coefficients from (29) are

[TABLE]

Figure 3 illustrates the accuracy of fulfillment of the boundary conditions (23)–(26). The absolute value of the difference between the exact solution (27) and the obtained approximate solution in the domain $D(s)$ is presented on Figure 4.

Note that for the problem considered one may look for an approximate solution in the form

[TABLE]

i.e., having only even coefficients $a_{2n}$. In such way the boundary condition (24) is satisfied automatically. The proposed algorithm can be applied with minimal modifications for such simplified form of an approximate solution, however we did not observe any gain in the obtained result, the original formulation of the algorithm performed equally well.

6 Possible extensions of the method to more general FBPs and final

remarks

6.1 Generalization of the operator $\mathbf{A}$

The presented method can be extended onto operators of the form

[TABLE]

with sufficiently regular coefficients $\alpha_{i}(x)$. For the construction of the formal powers in this case see [22]. The FBPs with this type of operators are very common in financial applications.

6.2 Generalizations of the conditions on a free boundary

6.2.1 Linear generalization

The conditions (7) and (8) are sometimes referred to as Stefan’s conditions. Our algorithm can be applied to more general conditions of the type

[TABLE]

for some first order linear differential operators $\mathbf{G}_{3}$ and $\mathbf{G}_{4},$ and where $g_{3}$ and $g_{4}$ are given functions of three variables.

6.2.2 Nonlinear generalization

For some FBPs we can only express the conditions (7) and (8) as

[TABLE]

where $g_{3}$ and $g_{4}$ are given functions of five variables. In this case the elegant decomposition into a linear and nonlinear problems will be lost, however the hierarchic structure of the minimization problem remains, i.e., the minimization Problem 8 can be reduced to the Problem 9 with the only difference that the auxiliary problem (16) may be nonlinear.

6.3 Nonlinear form of the boundary

Often in applications some additional information is available about the free boundary structure. With this additional knowledge (or for other reasons) the linear decomposition of the boundary might not be appropriate. It is easy to see that we can still apply the algorithm assuming that

[TABLE]

for the known function $\xi$. This nonlinear structure can slow down the algorithm for the search of the minimum in Problems 8 or 9, due to the constraints (15) and (17) respectively. In our particular Matlab implementation, we have used fmincon function that performs better with a linear constraint on the boundary than a nonlinear one.

6.4 Concluding remarks

A method for approximate solution of a large variety of FBPs is proposed. It is based on a possibility to construct a complete system of solutions of a parabolic equation called transmuted heat polynomials. The numerical implementation is relatively simple and direct. The time required for computations is within seconds. The method admits extensions onto a much larger class of FBPs then that discussed in the paper.

Acknowledgements

Research was supported by CONACYT, Mexico via the project 222478. The first named author would like to express his gratitude to the excellence scholarship granted by the Mexican Government via the Ministry of Foreign Affairs which gave him the opportunity to develop this work during his stay in the CINVESTAV, Mexico.