CrossMark

Available online at www.sciencedirect.com

ScienceDirect

Mathematics and Computers in Simulation 101 (2014) 61-77

Original Article

A meshless method for an inverse two-phase one-dimensional

nonlinear Stefan problem

B. Tomas Johansson a'b, Daniel Lesnicc, Thomas Reeved *

a EAS, Aston University, Birmingham B4 7ET, UK ' ITN, Campus Norrkoping, Linkoping University, Sweden c Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK d School of Mathematics, University of Birmingham, Birmingham B15 2TT, UK

Received 16 August 2012; received in revised form 17 May 2013; accepted 13 March 2014 Available online 29 March 2014

Abstract

We extend a meshless method of fundamental solutions recently proposed by the authors for the one-dimensional two-phase inverse linear Stefan problem, to the nonlinear case. In this latter situation the free surface is also considered unknown which is more realistic from the practical point of view. Building on the earlier work, the solution is approximated in each phase by a linear combination of fundamental solutions to the heat equation. The implementation and analysis are more complicated in the present situation since one needs to deal with a nonlinear minimization problem to identify the free surface. Furthermore, the inverse problem is ill-posed since small errors in the input measured data can cause large deviations in the desired solution. Therefore, regularization needs to be incorporated in the objective function which is minimized in order to obtain a stable solution. Numerical results are presented and discussed.

© 2014IMACS. Published by Elsevier B.V. All rights reserved. MSC: 35K05; 35A35; 65N35

Keywords: Heat conduction; Inverse Stefan problem; Method of fundamental solutions; Two-phase change; Regularization

1. Introduction

The two-phase direct Stefan problem requires determining the temperature distribution and the moving free interface when the initial and boundary conditions, as well as the thermal properties of the bi-material involved, are known, see e.g. [19]. Under various boundary conditions it was shown to be well-posed, see [3,5,6]. In contrast to the direct problem, inverse Stefan problems require determining some initial temperature and/or boundary conditions, and/or thermal properties from additional information, which may involve the partial knowledge of the free surface, the temperature measured at some points inside the medium or on the boundary, the heat flux, etc., see [10].

In comparison to the studies on the one-phase flow, the literature on solving two-phase inverse Stefan problems is much more scarce, see [1,17,20]. However, these inverse design Stefan problems were linear because in their

* Corresponding author. Tel.: +44 (0) 121 204 3978. E-mail addresses: tomas.johansson@liu.se (B.T. Johansson), reevet@maths.bham.ac.uk (T. Reeve).

http://dx.doi.org/10.10167j.matcom.2014.03.004

0378-4754/© 2014 IMACS. Published by Elsevier B.V. All rights reserved.

formulation the position of the moving interface is considered known. When the position of the moving interface is unknown and also no temperature or heat flux boundary conditions are specified on a part of the boundary then one deals with a nonlinear and ill-posed inverse Stefan problem, see [11]. Although the uniqueness of a solution in Holder spaces for such a class of two-phase inverse nonlinear Stefan problems holds, see [10,12], these problems are still ill-posed because there is no continuous dependence of the solution on the input data.

The plan of the paper is as follows. In Section 2, we give the mathematical formulation of the one-dimensional two-phase inverse nonlinear Stefan problem and point out its ill-posedness. In Section 3, we describe the regularized numerical meshless method of fundamental solutions (MFS) for constructing a stable solution to the inverse problem. Section 4 presents and discusses numerical results obtained for some typical test examples with and without noise included in the input data. Finally, Section 5 presents conclusions and possible future work.

2. Mathematical formulation

Assume that the interface

s(t) e (0, I], for t e (0, T], (2.1)

and i(0) is given, and denote by

Qt = {(x, t) e (0,1] x (0, T]} the two-phase rectangular domain (0,1) x (0, T], which is subdivided by the interface into the two subdomains D\ = {(x, t) e QT|0 < x < s(t), t e (0, T]},

DT = {(x, t) e Qt|s(t) < x < I, t e (0, T]}.

We investigate the inverse nonlinear two-phase one-dimensional Stefan problem which requires finding the triplet solution (u1,u2,s) e C2,1(D1t) x C^D) x (C([0, T]) nC'(0, T])), satisfying (2.1), the heat equations

du i d2ui

~dT = 01 Ix2 '

(x, t) e DT (2.2a)

duo d2U2 -,

Ht =a2, {xJ) e °T (2'2b)

where an > 0 is the thermal diffusivity of the heat conductor DT for n =1,2, the initial conditions at t = 0

ul(x, 0) = u0(x), x e [0, i(0)], (2.3a)

U2(x, 0) = u2(x), x e [s(0), l], (2.3b) the interface Stefan conditions at x=s(t)

u1(s(t),t) = u2(s(t),t) = u*(t), t e (0, T], (2.4a) , dui du2

s'(t) =-Ki --1(s(t),t) + K2--2(s(t),t), t e (0, T], (2.4b) dx dx

where Kn = kn/(p2L), and kn, pn, L are the thermal conductivities, densities and latent heat, respectively, of DlT (water) and DT (ice), and the Cauchy boundary conditions at x=l

u2(l, t) = fi (t), t e (0, T], (2.5a)

k2dui(l, t) = gi(t), t e (0, T], (2.5b) dx

We also impose the compatibility conditions at the corners (x, t) = (s(0), 0) and (l, 0), namely

u1(s(0)) = u0(s(0)) = u*(0), (2.6a) r\ du 0

u0(l) = fi(0), k2 ^(0) = gi(0). (2.6b)

Data (2.7a) and (2.7b) missing

(2.4a) !

(2.4b)./ N 2 v .> (2.2b) in D|, (2.5a)

(2.2a) in DT ls(t) = ? (2.5b)

; (unknown)

(2.3a) s(0) (2.3b) x

Fig. 1. Sketch of the two-phase inverse nonlinear Stefan problem (2.1)—(2.5).

Remark that the fixed boundary x = l is overspecified because both Dirichlet and Neumann (i.e. Cauchy) data are prescribed in (2.5a) and (2.5b), whilst the fixed boundary x=0 is underspecified because no condition is prescribed on it. Physical quantities of interest related to this hostile or inaccessible boundary are the temperature

ui(0,t) = fo(t), t e (0, T], (2.7a)

the heat flux dui

-ki -^(0, t) = go(t), t e (0, T], (2.7b)

and the mass

r s(t)

/ u1(x, t) dx = E(t), t e (0, T], (2.7c)

This data can also be required to satisfy the compatibility conditions at the origin (x, t) = (0, 0), namely,

0 dui fs(0) 0

f0(0) = u?(0), g0(0) = -ki -±(0), E(0) = u0(x) dx, (2.8)

A sketch of the conditions in the two-phase inverse nonlinear Stefan problem is shown in Fig. 1. Sometimes, in practice, the measurement of both the boundary temperature (2.5a) and the heat flux (2.5b) may not be easy and only one of this data may be available. In such a situation one could measure instead the upper base (i.e. final time) internal temperature at t = T, namely

ui(x, T) = u\(x), x e [0, s(T)], (2.9a)

u2(x, T) = uT(x), x e [s(T), I], (2.9b)

where s(T) e (0, l) is also given. However, it turns out, see [11], that this latter 'upper base data (2.9)' inverse problem is more ill-posed than the former 'Cauchy data (2.5a and 2.5b)' inverse problem because it can have more than one solution, i.e. the solution is not unique. However, as we shall see below, the former problem can have at most one solution, i.e. the solution is unique, though the problem is still ill-posed since small errors in the input data (2.5) can cause large errors in the output data (2.7). Therefore, in this study only the inverse problem given by Eqs. (2.1)-(2.6) will be investigated. This problem may be considered as a continuation problem of the solution of the parabolic linear heat equation from the boundary x = l, where the Cauchy data (2.5a and 2.5b) is given, into the domain QT. It can also be reinterpreted as a "backward in space" inverse heat conduction problem with an "initial" transient boundary, see [18]. However, in contrast to the non-characteristic Cauchy problem, there is the unknown phase transition moving boundary x = s(t) in QT also to be determined and this essentially complicates the task of analytic continuation.

An initial attempt would be to split the two-phase inverse Stefan problem (2.1)-(2.6) into two problems. The first is the nonlinear inverse boundary determination problem for the pair (u2, s) satisfying equations (2.1), (2.2b), (2.3b), (2.5a and 2.5b), (2.6) and

u2(s(t), t) = u*(t), t e (0, T]

(2.10)

M°(s(0)) = u*(0). (2.11)

The solution of this problem is unique, see [9], even if the initial condition (2.3b) is not imposed, see [22]. Once the boundary x=s(t) and the heat flux dur(s(t), t) have been determined, the second is the linear inverse problem for determining the temperature u1 satisfying equations (2.2a), (2.3a),

u1(s(t), t) = u*(t), t e (0, T] (2.12)

u01(s(0)) = u*(0), (2.13)

dui , du2

-Ki--1(s(t),t) = s'(t) -K2-T2(s(t), t), t e(0,T]. (2.14)

The solution of this problem is unique, see [4] even if the initial condition (2.3a) is not imposed, see [14,21]. However, both the above problems are ill-posed since their solutions do not depend continuously on the input data. In order to obtain stable solutions regularization is necessary and this involves choosing at least two regularization parameters, one for each problem. Besides depending on the amount of noise in the input data they depend on each other since the first nonlinear ill-posed problem has to be solved first to provide the input for the second linear ill-posed problem. So, this two-parameter choice becomes complicated.

Therefore, it appears more useful to solve the composite problem (2.1)-(2.6) in one go for simultaneously determining the solution (u1, u2, s) .

The problem is still ill-posed but it can now be regularized by choosing a single regularization parameter for the unknown temperature fields. This combined approach has been used previously by the authors for simultaneously determining a heat source and the initial temperature, see [15]. Finally, we will also investigate the case when the initial conditions (2.3a) and (2.3b), together with the compatibility conditions (2.6a) and (2.6b), are not prescribed.

3. The method of fundamental solutions

We approximate the solutions (u\, u2) of the heat equations (2.2a) and (2.2b) using the method of fundamental solutions (MFS) for the unsteady heat conduction in composite layered materials, recently developed by the authors in [7,16,17]. In the MFS, we approximate ui and u2 by linear combinations of fundamental solutions of the heat equation

Gn(x, t; y, t) =

H(t — t) ^J4ann(t — t )

(x — y)2 4an(t — t)

n =1, 2,

where H is the Heaviside function, of the form M

(x, t) = c(pGn(x, t; y{f\ Tj), (x, t) eDnT, n = 1, 2.

In expression (3.2), the source points (y(n), Tj) for j = 1, M and n = 1, 2, are located outside the solution domains Dj-in the following way:

j = 1, 2M1

yf =\ h + S3M1—j+1, j = 2M1 + 1, 3M1

h + S4M1 —j+1, j = 3Mx + 1, 4M1

—h + SM1—j+1, j =1,M1

y? ={—h + S2M1 —j+1, j = M1 + 1, 2M1

(3.3a)

(3.3b)

I + h,

j = 2Mi + 1, 4Mi

Fig. 2. Representation of the two-phase problem with domains DT and D2T, unknown boundary data (x), and source points (•••) and (---) placed

externally to the domains DT and D^, respectively.

where h is a preassigned MFS positive parameter giving the distance between the source points and the boundary, M = 4M1,

(2j - 1 - 2M0T

for j = 1, 2Mi,

Tj = Tj-2Ml for j = 2M1 + 1,4M1; and sj = s(j+M1 ) for j = 1, M1, see Fig. 2 for a representation of the domains, boundaries and placement of source points.

(n) n=1,2

The2M unknown coefficients c = (cj ' and the M1 = M/4 discretization points s = (sj )j=1M ofthe unknown free boundary have to be determined by collocating the conditions (2.3)-(2.6), as described next. Let us select a uniform distribution of collocation points given by

(2i — 1)T _

t0 = 0, ti =- = |T2M1-i+1|, i =1,M1,

~ti = — for i = 0, M2,

' M2 , 2

(k) ks(0) (k) ... *(/ — s(0)) x\ = —--, x2 = s(0) +

k =1,K,

and collocate equations (2.3)—(2.6) as follows:

Un(x(k), 0) = u0n(x(nky), k = 17K, n =1, 2,

U1(si, ti) = U2(si, ti) = u (ti), i = 0, M1,

s(t1) — s(0)

du1 du2 ,

—K1 d^(s!,t!) + K2 --2(sl,ti) = s'(ti) = )

dx dx I s(ti) — s(ti—1)

ti — ti 1

du 2 ,

U2(l, ~ti) = fi(ti), k2 — (l, ~ti) = gi(ti), i = 0, M2. dx

i = 1 , M1 ,

i = 2, M1

In total, via (3.2), Eqs. (3.5)-(3.8) form a system of (2K +3M1 +2M2 + 4) equations with 9M1 unknowns (c, s). In general, we require to have at least as many equations as unknowns and therefore we require

2M2 + 2K + 4>6M1, or M2 + K + 2>3M1.

Note that this system is linear in c, but it is nonlinear in s. In addition, our inverse problem is ill-posed. Therefore, we apply the nonlinear Tikhonov regularization method based on minimizing the non-linear regularized least-squares functional

2 K M \ Mi M

F(c, s) = ]T]T £4n)G„(*® 0; yjn), Tj ) - u^®) + E E^ife j ) - u\t¡)

n=i k=i\ j=i ) i=0\ j=i

Mi ! m

+ E Y,c?G2(si, ti; y?, Tj) - u*(ti) i=0\ j=i

i=i V j=i

Kj ^ t; yf, j ) - Kij 'G-iSi, tt; j, t; )

- s'(ti)

+E I EC/2)G2(1, tt; j, Tj ) - Mi)] Ecj2)k2 ^(t u; j ) - gm

i=0 \j=i

i=0 j=i

+ Ài|c|2 + À2|s|2,

where A.i, >0 are regularization parameters which can be prescribed according to some criterion, e.g. the L-surface, or, more simply, by trial and error. Note that in the last term of (3.9), for simplicity and in a first attempt, we have imposed that s e C[0, T], but more regularity such as s e C1[0, T] can also be imposed in order to get some stability estimates.

In imposing the flux boundary conditions (2.4b) and (2.5b) the x-partial derivative of (3.1) is needed, as given by dGn (x - y)H(t - t) ( (x - y)2 \

—(x,t;y,T) =----exH ^^-: h n = 1 2-

dx 2^4na3n(t - t)3 v 4a"(t -Tv

The minimization of the functional (3.9) is performed using the MATLAB toolbox Isqnonlin which is designed to minimize a sum of squares of arbitrary differentiable functions. The gradient does not need to be supplied by the user and the simple bounds on the variables

0 <si<l for i =i,Mi (3.i0)

are also allowed.

The initial guess to start the iterative process is arbitrary and in this study we take c0 = 0 and s0 = s(0)

4. Numerical results and discussion

In this section, we apply the MFS outlined in the previous section to the inverse two-phase nonlinear Stefan problem (2.1)-(2.6). In the first two examples we have analytical solutions available and we also investigate for one of them the case when the initial conditions (2.3a) and (2.3b) are not given. In the third example, an analytical solution is not available. Previous experience with applying the MFS for the heat equation, [16,17], suggests that the MFS parameter h should be not too small (which will result in less accurate approximations) nor too large (which will increase the ill-conditioning of the system). In this section, the value of h, as well as the values of the regularization parameters and k2, are chosen by trial and error. Nevertheless, more rigorous investigations on these choices should be undertaken in a future work.

4.1. Example 1

Let us take T =1, l = 2, «1 = 2, k1 = 1, a2 = 1, fe = 2, ^1 = 1, ^2 = 2. We also choose in Eqs. (2.3a) and (2.3b) the initial conditions

u1(x, 0) = u1(x) = 2

0.5— x exp —Ô— — 1

x e [0, s(0) = 0.5] (4.1a)

u2(x, 0) = u2(x) = exp(0,5 - x) - 1, x e [0,5 = s(0), I = 2], (4.1b) We also take u*(t) = 0 such that (2.4a) becomes

u1(s(t),t) = u2(s(t),t) = 0, t e (0, T = 1], (4.2) and take the Cauchy boundary conditions (2.5a) and (2.5b) specified by

u2(2,t) = f(t) = exp(t -1,5) - 1, t e [0, T = 1], (4.3a)

k2du2(2,t) = gi(t) = -2exp(t -1,5), t e [0, T = 1], (4.3b)

Then the analytical solution of the inverse Stefan problem (2.1), (2.2), (2.4), (4.1)-(4.3) is given by

u1(x,t) = 2^exp^t + °25 -x^j - , (x,t) e ~D\, (4.4a)

u2(x,t) = exp(t + 0,5 -x) - 1, (x,t) e DT, (4.4b)

s(t) = t + 0,5, t e [0, T = 1], (4.4c)

which can be verified by direct substitution.

Of particular interest is to recover the missing data at the inaccessible hostile boundary x=0, given by the boundary temperature

u1(0,t) = f0(t) = 2^exp(, t e[0, T = 1] (4.5a)

the heat flux

du 1 (1 + 0,5 N

-k1 -¿(0't) = g0(t) = expi —J, t e [0, T = 1] (4.5b)

and the mass

[s(t) ( (t + 0,5 \ N

J u1(x,t) dx = E(t) = 2i2expi—-— j-t - 2,5 J , t e [0, T = 1], (4.5c)

If we consider the boundary temperature (4.3a) at x=l=2 as the quantity which is physically measured in practice, whilst the heat flux (4.3b) is prescribed, in order to take into account the uncertainties in the measurement we add to it random additive noise generated as

u2(1, t) = f(t) = f (t) + N (0, a2), (4.6)

where N(0, a2) represents the normal distribution with mean zero and standard deviation a which is taken to be

a = 5 x max |f(t)|, (4.7)

te[0,1]

where 5 is the percentage of noise. In MATLAB, N(0, a2) is generated by multiplying a with randn (a matrix of normally distributed pseudorandom numbers).

Fig. 3(a)-(d) shows the numerical MFS solutions for s(t), u1(0, t), -k1 d|L(0, t) and E(t), as functions of t e [0, 1], obtained with M1 = 20, M2 = 40, K = 40, i.e. 2K + 3M1 + 2M2 + 4 =224 equations with 9M1 = 180 unknowns, for no

0.75 1

1.25 1.5

0 0.2 0.4 0.6 0.8 1

Fig. 3. The numerical (•••) MFS solutions for: (a) s(t), (b) u1(0, t), (c) —k1 duL(0, t), and (d) E(t), as a function of t e [0, 1], obtained with M1 = 20, M2 = 40, K = 40, h = 3 for no noise S = 0 and À1 = À2 = 0 in comparison with the exact (—) solutions (4.4c), (4.5a)-(4.5c), respectively, for Example 1.

noise, i.e. 5 = 0, and = k2 = 0, in comparison with the exact solutions (4.4c), (4.5a)-(4.5c), respectively. In calculating the mass (2.7c) using the MFS the following integral is needed:

, ^ j H (t - t ) I" f( s(t) -y \ f( -y

G1(x,t; y,T) dx =- erf , = - erf , = , (4.8)

J0 ' 2 L VV4«1(t t)) VV4«1(t -r)J\ '

where erf is the error function.

Fig. 4(a)-(d) shows the same as Fig. 3(a)-(d), but for 5 = 1% and 5 = 5% noise and = k2 = 10-6. In Figs. 3 and 4 it can be seen that stable and accurate numerical MFS solutions are obtained.

Next, we look at the case when the initial conditions (4.1a) and (4.1b) are not imposed. In this case the first term in the right-hand side of (3.9) drops out, i.e. Eq. (3.5) is not imposed, such that, via (3.2), Eqs. (3.6)-(3.8) form now a smaller system of (3M1 + 2M2 + 4) equations with 9M1 unknowns, and we require M2 + 2 >3M1.

Figs. 5(a)-(d) and 6(a)-(d) represent the same quantities as Figs. 3(a)-(d) and 4(a)-(d), respectively, but obtained with M1 = 10 and M2 = 40, i.e. 3M1 + 2M2 + 4 = 114 equations with 9M1 = 90 unknowns. In addition, Figs. 5e and 6e present the numerical results for the initial temperature

f u 1 (x, 0), x e [0, s(0) = 0,5] u(x, 0) = I (4.9)

[ u2(x, 0), x e [0,5 = s(0), I = 2]

in comparison with the exact solution given by Eqs. (4.1a) and (4.1b). Again, even when we are missing the initial data (2.3a) and (2.3b), we obtain accurate and stable results in all figures, although some instabilities start to manifest in the retrieved heat flux in Fig. 6(c).

So far, Example 1 has investigated retrieving a linear function of t for the free surface (4.4c). In the next example we consider a nonlinear variation.

Fig. 4. The numerical MFS solutions for: (a) s(t), (b) «1(0, t), (c) -k1 ^(0, t), and (d) E(t), as a function of t e [0, 1], obtained with M1 = 10, M2 =20, K=20, h = 3 for S = 1% (•) and S = 5% (■) noise and X1 = X2 = 10-6 in comparison with the exact (—) solutions (4.4c), (4.5a)-(4.5c), respectively, for Example 1.

4.2. Example 2

Let us take T = 1, l = 1, a\ =0.5, a2 = 1, ki =0.5, k2 = 1, L = 1, p2 = 1, K =0.5, K2 = 1, and u* =0. We also choose in Eqs. (2.3a) and (2.3b) the initial conditions

erf erf

0 v 2Vait0 ,

ui(x, 0) = u01(x) =1--\v 1 , xe[0,s(0) = 0.4] (4.10a)

erfc erfc

0 V 2Ja2t0 ,

u2(x, 0) = u0(x) =1 +-v,v 2 , x e [0.4 = s(0),l = 1] (4.10b)

where y = 0.479611, t0 = 0.695571 and erf and erfc are the error and the complementary error function, respectively. We also take the interface condition (4.2) and the Cauchy boundary conditions (2.5a) and (2.5b) specified by

erfc ' 1

. 2 Ja2(t+t0) .

u2(1,t) = fi(t) = -1 +-v ^ 2 , V , t e (0, T =1], (4.11a)

erfc №

k2 exp (-40^)

du2 , 4a,(t+t0) ,

k2-z1(1,t) = gi(t) =--, v 2( +T^Y, t e (0, T =1], (4.11b)

dx vWt + t0) erfc (^J

0.4 0.6

Fig. 5. The numerical (•••) MFS solutions for: (a) s(t), (b) u1(0, t), (c) -k1 ^7(0, t), and (d) E(t), as a function of t e [0, 1], and (e) u(x, 0) as a function of x e [0, 2], obtained with M1 = 10, M2 = 40, h = 3 for no noise 5 = 0 and = ^2 = 0 in comparison with the exact (—) solutions (4.4c), (4.5a)-(4.5c), respectively, for Example 1.

Then the analytic solution of the inverse Stefan problem (2.1), (2.2), (2.4), (4.2), (4.10) and (4.11) is given by, see [7,8,17],

u\(x, t) =1 —

2 V a1(t+t0)

(x, t) e DT,

(4.12a)

U2(x, t) = — 1 +

2 V a2(t+t0)

(x, t) e DT,

s(t) = y\/t + t0, t e [0, T = 1],

(4.12b)

(4.12c)

Fig. 6. The numerical MFS solutions for: (a) s(t), (b) u1 (0, t), (c) -k1 dU1(0, t), and (d) E(t), as a function of t e [0, 1], and (e) u(x, 0) as a function of x e [0, 2], obtained with M1 = 10, M2 = 40, h = 3 for S = 1% (•) and S = 5% (■) noise with = A.2 = 10-6 in comparison with the exact (—) solutions (4.4c), (4.5a)-(4.5c), respectively, for Example 1.

which can be verified by direct substitution. Of particular interest is to recover the missing data at the inaccessible hostile boundary x = 0, given by the boundary temperature

u1(0,t) = f0(t) = 1, t e [0, T = 1] (4.13a)

the heat flux

-ki ~r1(0,t) = g0(t) = 1 -—, t e [0, T =1] (4.13b)

dx vWt + t0) erf^

and the mass

s(t) 2 /ax(t + t0) / ( y2 , ,

u1(x,t) dx = E(t) =---u ; 1 -exp - — , t e [0, T = 1]. (4.13c)

0 erf ( tt^M n V V 4a1

K2y/ar

Noise is added in the boundary temperature (4.11a), as in (4.6).

0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 t

0.2 0.4 0.6 0.8 1

S- 0.95 » 0.9 0.85 0.8

0.3 0.28

0.24 0.22 0. 0.18

0.4 0.6

0.4 0.6

Fig. 7. The numerical (•••) MFS solutions for: (a) s(t), (b) u1(0, t), (c) —k1 ^(0, t), and (d) E(t), as a function of t e [0, 1], obtained with M1 =20, M2 =40, K=40, h = 2 for no noise S = 0 and À1 = À2 = 0 in comparison with the exact (—) solutions (4.12c), (4.13a)-(4.13c), respectively, for Example 2.

Fig. 8. The numerical MFS solutions for: (a) s(t), (b) u1(0, t), (c) —k1 duL(0, t), and (d) E(t), as a function of te [0, 1], obtained with M1 =20, M2 = 40, K=40, h = 2 for S = 1% (•) and S = 5% (■) noise and À1 = À2 = 10—6 in comparison with the exact (—) solutions (4.12c), (4.13a)-(4.13c), respectively, for Example 2.

0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 t

0.4 0.6

0.4 0.6

0.4 0.6

Fig. 9. The numerical (•••) MFS solutions for: (a) s(t), (b) ui(0, t), (c) -k1 dU1(0,t), and (d) E(t), as afunctionof t e [0, 1], obtained with M1 =40, : 80, h = 2 for no noise S = 0 and À1 = À2 = 10-12 in comparison

Example 2.

M2 = 80, K = 80, h = 2 for no noise S = 0 and X1 = X2 = 10 12 in comparison with the exact (—) solutions (4.12c), (4.13a)-(4.13c), respectively, for

Figs. 7 and 8 for Example 2 are the analogous of Figs. 3 and 4 for Example 1 and the same conclusions about the good performance of the method are obtained. To further improve the results obtained in Fig. 7, in Fig. 9 we increase the number of collocation and source points.

In both Examples 1 and 2 analytic solutions were available and this enabled assessing the accuracy of the numerically obtained results. The next example considers the case when an analytic solution is not available.

4.3. Example 3

In this example, partially taken from [1] where a linear inverse Stefan problem was addressed, an analytic solution is not available. We take T = 1, l = n/2, «1=2, «2=1, s(0) = n/4, ^1 = 1, K2 = 2, k1 = 1, k2 = 2, u* = 0, gl = 0, and

Fig.10. The direct problem MFS for u2 (n/2, t) (—), as a function of t e [0,1], obtained with h = 2, M1 =40, M2 =60, K=58 (i.e. K+M1 + M2 + 2=160 equations with 4M1 = 160 unknowns) and X = 10-14, and the perturbed data for S =1% (■) and S = 5% (A), for Example 3.

u0(x) = cos(2x) for x e [|, §]. Since an analytic solution is not available the data (2.5a) is numerically simulated by solving separately using the MFS the direct problem in the domain D\ given by Eqs. (2.2b), namely

du2 d2u2

HT = d2,

(x, t) e DT = (s(t), l) X (0, T],

Eq. (2.3b) given by

u2(x, 0) = u2(x) = cos(2x), x e condition (2.4a) given by

u2(s(t),t) = u*(t) = 0, t e (0, 1],

Û , 2 J

(4.14)

(4.15)

(4.16)

0.4 0.6 0.8

S 1.5 1

Fig. 11. The numerical MFS solutions (•••) for: (a) s(t), (b) u1 (0, t), (c) -k1 ddLL(0, t), and (d) E(t), as a function of te [0, 1], and (e) u1 (x, 0), as a function of x e [0, n/4] obtained with M1 = 20, M2 = 40, K = 40, for no noise 5 = 0 and = ^2 = 0, for Example 3. In case (a), the comparison with the exact solution (4.18) for s(t) is also shown, and in cases (b)-(e) comparisons are made with the MFS approximation (—) that was generated in [17].

the boundary condition (2.5b) given by du2 tn \

k2(2, V = g'(t) = 0, t e (0, 1],

when the free surface is known and is given by

s(t) = arctan(1 + t), t e (0, 1], Observe that the compatibility conditions in (2.6a) and (2.6b) given by

u2(s(0)) = u*(0), k2d2 = gi(0)

(4.17)

(4.18)

(4.19)

-—) 2.5

Si 8 2

4 3.5 3

,—) 2.5

^ 1.5 1

0 0.2 0.4 0.6

0.4 0.6 0.8

Fig. 12. The numerical MFS solutions for: (a) s(t), (b) u1 (0, t), (c) -k1 u (0, t), and (d) E(t), as a function of t e [0, 1], and (e) u1 (x, 0), as a function of x e [0, n/4] obtained with M1 =20, M2 =40, K=40, for 5 = 1% (■) noise and l1 = X2 = 10-6, for Example 3. In case (a), the comparison with the exact solution (4.18) for s(t) is also shown, and in cases (b)-(e) comparisons are made with the MFS approximation (—) that was generated in [17].

are automatically satisfied by the data (4.15)-(4.18). In the direct problem we collocate the Eqs. (4.15)-(4.17) and (4.19) as

u2(x2), 0) = w0(x2k)), k = T7K (4.20)

U2(si, ti) = u*(ti), i = 0, M\ (4.21)

k2~±(l, ti) = gi(ii), i = 0, M2 (4.22)

resulting in, via (3.2), a system of (K+Mi + M2 + 2) linear equations with 4Mi unknowns c(2) = (cj2)) .^-m. A nec-

J J—1,M

essary condition for a unique solution is K+M2 + 2 >3M1. This system of equations is ill-conditioned and therefore we employ linear Tikhonov regularization with regularization parameter X. The numerical results obtained for the boundary temperature u2(n/2, t), as a function of t e [0, 1], are shown in Fig. 10. This data, to which noise is added as in (4.6), is then used as the input (2.5a) in the inverse problem.

The imposition of the initial condition (2.3a) cannot be arbitrary and in order to ensure that the solution exists we simply do not impose it. In this case the first term in the right-hand side of (3.9) for n = 1 drops out, i.e. Eq. (3.5) for n = 1 is not imposed, such that, via (3.2), equations (3.5) for n = 2, (3.6)-(3.8) form a system of (3M1+ 2M2 + K+4) equations with 9M1 unknowns, and we require 2M2 + K+4 >6M1.

Figs. 11(a)-(e) and 12(a)-(e) show the MFS solutions for s(t), u1(0, t), -k1 ^^(0, t), E(t) and u1(x, 0) obtained for M1 = 20, M2 =40, K =40 (i.e. 3M1 + 2M2 + K+4 =184 equations with 9M1 = 180 unknowns) for S = 0 and X1 = X2 = 0, S = 1% and X1 = X2 = 10-6. We note that to generate the stable results in Fig. 11 the variable and function tolerances used in the MATLAB toolbox Isqnonlin were reduced from 10-6 to 10-9. For the free boundary, comparison between the MFS solution and the exact solution (4.18) made in Fig. 11(a) shows good agreement. We point out that the reconstructions become rather inaccurate when increasing the noise level, and for S = 5% it is not possible to obtain a reasonable approximation. Moreover, by inspecting Fig. 12(a)-(e) (see also Figs. 6(c) and 8(c)) we note that the numerical solutions seem to deviate from the 'exact' MFS approximations that were generated in [17], as the time t approaches the final time T. This is to be somewhat expected since if we want to recover the heat flux at x=0 over the whole time interval [0, T] then we need to use the Cauchy data (2.5a) and (2.5b) over an extended interval [0, T + r], where r >0 is related to the concept of 'future times' in the inverse heat conduction literature, see [2,13].

5. Conclusions

In this paper, the one-dimensional two-phase nonlinear inverse Stefan problem has been investigated using a regularized MFS. For both exact and noisy data the method has been shown to be accurate, stable and robust. At present there are no other results available to compare our MFS with since this paper undertakes the first numerical investigation to solve the inverse two-phase nonlinear Stefan problem. We mention that alternatively, one could employ the boundary element method, as a powerful and well-suited numerical boundary discretisation approach for solving the problem. Future work will concern extending the MFS developed in this study to multi-dimensional nonlinear Stefan problems.

Acknowledgements

T. Reeve would like to thank the financial support received from the EPSRC. B.T. Johansson and D. Lesnic would like to thank the Isaac Newton Institute for the Visiting Fellowships within the programme on "Inverse Problems" offered in November-December 2011 during which period some of this work has been performed. The comments and suggestions made by the referees are acknowledged.

References

[1] D.D. Ang, A. Pham Ngoc Dinh, D.N. Tranh, Regularization of an inverse two-phase Stefan problem, Nonlinear Anal. 34 (1998) 719-731.

[2] J.V. Beck, B. Blackwell, C.R. St-Clair Jr., Inverse Heat Conduction, J. Wiley-Intersc. Publ., New York, 1985.

[3] J.R. Cannon, D.B. Henry, D.B. Kotlow, Classical solutions of the one-dimensional, two-phase Stefan problem, Ann. Mat. Pura Appl. 107 (1975)311-341.

[4] J.R. Cannon, J. Douglas Jr., The Cauchy problem for the heat equation, SIAM J. Numer. Anal. 4 (3) (1967) 317-336.

[5] J.R. Cannon, M. Primicerio, A two phase Stefan problem with flux boundary conditions, Ann. Mat. Pura Appl. 88 (1971) 193-205.

[6] J.R. Cannon, M. Primicerio, A two phase Stefan problem with temperature boundary conditions, Ann. Mat. Pura Appl. 88 (1971) 177-191.

[7] S. Chantasiriwan, B.T. Johansson, D. Lesnic, The method of fundamental solutions for free surface Stefan problems, Eng. Anal. Bound. Elem. 33 (2009) 529-538.

[8] J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, 1984.

[9] A. El Badia, F. Moutazaim, A one-phase inverse Stefan problem, Inverse Prob. 15 (1999) 1507-1522.

[10] N.L. Gol'dman, Inverse Stefan Problems, Kluwer Academic Publ., Dordrecht, 1997.

[11] N.L. Gol'dman, Inverse problems with phase transitions., in: O.M. Alifanov, N.V. Kerov, V.V.Michailov, A.N. Nenarokomov (Eds.), Proc. 3rd Int. Conf. Dynamic Systems Identification and Inverse Problems, Moscow-St. Petersburg, Russia, 1998, pp. 128-139.

[12] N.L. Gol'dman, Properties of solutions of the inverse Stefan problem, Differ. Equ. 39 (2003) 66-72.

[13] D. Nho Hao, Methods for Inverse Heat Conduction Problems, Peter Lang, Frankfurt am Main, 1998.

[14] C.D. Hill, Parabolic equations in one space variable and the non-characteristic Cauchy problem, Comm. Pure Appl. Math. 20 (1967) 619-633.

[15] B.T. Johansson, D. Lesnic, A procedure for determining a spacewise dependent heat source and the initial temperature, Appl. Anal. 87 (2008) 265-276.

[16] B.T. Johansson, D. Lesnic, A method of fundamental solutions for transient heat conduction in layered materials, Eng. Anal. Bound. Elem. 33 (2009) 1362-1367.

[17] B.T. Johansson, D. Lesnic, T. Reeve, A meshless method for an inverse two-phase one-dimensional linear Stefan problem, Inverse Prob. Sci. Eng. 21 (1) (2013) 17-33.

[18] D.A. Murio, Solution of inverse heat conduction problems with phase changes by the mollification method, Comput. Math. Appl. 24 (1992) 45-57.

[19] L. Rubinstein, The Stefan problem: Comments on its present state, J. Inst. Maths Appl. 24 (1979) 259-277.

[20] D. Slota, Solving the inverse Stefan design problem using genetic algorithms, Inverse Prob. Sci. Eng. 16 (2008) 829-846.

[21] Y.B. Wang, J. Cheng, J. Nakagawa, M. Yamamoto, A numerical method for solving the inverse heat conduction problem without the initial value, Inverse Prob. Sci. Eng. 18 (2010) 655-671.

[22] T. Wei, M. Yamamoto, Reconstruction of a moving boundary from Cauchy data in one-dimensional heat equation, Inverse Prob. Sci. Eng. 17 (2009) 551-567.