# The Solution of Two-Phase Inverse Stefan Problem Based on a Hybrid Method with Optimization

Authors: Yang Yu; Xiaochuan Luo; Haijuan Cui
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The two-phase Stefan problem is widely used in industrial field. This paper focuses on solving the two-phase inverse Stefan problem when the interface moving is unknown, which is more realistic from the practical point of view. With the help of optimization method, the paper presents a hybrid method which combines the homotopy perturbation method with the improved Adomian decomposition method to solve this problem. Simulation experiment demonstrates the validity of this method. Optimization method plays a very important role in this paper, so we propose a modified spectral DY conjugate gradient method. And the convergence of this method is given. Simulation experiment illustrates the effectiveness of this modified spectral DY conjugate gradient method.


---

## Body

## 1. Introduction

The process of solidification/melting occurs in many industrial fields, such as steel-making continuous casting. In the continuous casting process [1, 2], molten steel comes into slab by solidification. And this process is shown in Figure 1. Because of the heat symmetry, we only consider 1/2 of slab thickness, and we divide this domain into two regions by the freezing front, which is shown in Figure 2. D
1 stands for the solid phase and D
2 is taken by the liquid phase. The solidification problem of continuous casting is associated with Stefan problem, which is a model of the process with a phase change.Figure 1
The solidification process of continuous casting slab (x stands for the direction of heat transfer and z stands for casting direction).Figure 2
The domain formulation of the two-phase Stefan problem.Many researches have focused on solving Stefan problems [3–7]. Grzymkowski and Slota [8, 9] applied the Adomian decomposition method (ADM) combined with optimization for solving the direct and the inverse one-phase Stefan problem, and Hetmaniok et al. [10] used the Adomian decomposition method (ADM) and variational iteration method to solve the one-phase Stefan moving boundary problems. In comparison to the studies on the one-phase Stefan problems, researches on two-phase Stefan problems [11] are much more scarce. Słota [12] applied homotopy perturbation method (HPM) to solve the two-phase inverse Stefan problem. However, in their formulation, these two-phase inverse Stefan problems were solved when the position of the moving interface is known. While the position of the moving interface is unknown and no temperature or heat flux boundary conditions are specified on one part of the boundary, this two-phase inverse Stefan problem may be difficult to be solved. This problem is close to the problem which occurs in the continuous casting. Therefore, this paper focuses on solving this problem by a hybrid method with optimization which combines the homotopy perturbation method (HPM) with the improved Adomian decomposition method (IADM).This paper is organized as follows. In Section2, we give the description of the two-phase Stefan problem. And the hybrid method for solving this problem is introduced in Section 3. We give a modified spectral DY conjugate gradient method and the convergence of this method in Section 4. In Section 5, some numerical results are given to illustrate the validity of the hybrid method.


## 2. The Description of Two-Phase Inverse Stefan Problem

In this section, the two-phase inverse Stefan problem is considered. The two domains are described (see Figure2):(1)
D
1
=
x
,
t
;
x
∈
0
,
ξ
t
,
t
∈
0
,
t
∗
,
D
2
=
x
,
t
;
x
∈
ξ
t
,
d
,
t
∈
0
,
t
∗
,and their boundaries with unknown function x
=
ξ
(
t
) are(2)
Γ
1
=
x
,
0
;
x
∈
0
,
s
,
s
=
ξ
0
,
(3)
Γ
2
=
x
,
0
;
x
∈
s
,
d
,
s
=
ξ
0
,
(4)
Γ
3
=
0
,
t
;
t
∈
0
,
t
∗
,
(5)
Γ
4
=
d
,
t
;
t
∈
0
,
t
∗
,
(6)
Γ
5
=
x
,
t
;
t
∈
0
,
t
∗
,
x
=
ξ
t
,
(7)
Γ
6
=
x
,
t
∗
;
x
∈
0
,
ξ
t
∗
.In practice, the measurement of temperature on the boundaryΓ
4 may not be easily obtained and the position of the moving interface described by means of function x
=
ξ
(
t
) is unknown. Therefore, the designed two-phase inverse Stefan problem consists in the calculation of the temperature distribution in the domains D
1 and D
2, respectively, as well as in the reconstruction of the functions describing moving interface x
=
ξ
(
t
). In order to solve this inverse problem, the additional information is given. In such a situation, the final time internal temperature at t
=
t
∗ on the boundary Γ
6 is known, which can be measured in practice. The temperature distribution is determined by the functions u
1
(
x
,
t
) and u
2
(
x
,
t
), which are defined in the domains D
1 and D
2, respectively, and they satisfy the following heat conduction equations:(8)
∂
u
1
x
,
t
∂
t
=
α
1
∂
2
u
1
x
,
t
∂
x
2
in
D
1
,
(9)
∂
u
2
x
,
t
∂
t
=
α
2
∂
2
u
2
x
,
t
∂
x
2
in
D
2
,where t and x refer to time and spatial location and α
1 and α
2 are the thermal diffusivity in the solid phase and the liquid phase, respectively. On the boundaries Γ
1 and Γ
2, the initial conditions are described in the following:(10)
u
1
x
,
0
=
φ
1
x
on
⁡
Γ
1
,
(11)
u
2
x
,
0
=
φ
2
x
on
⁡
Γ
2
,
(12)
ξ
0
=
s
.On the boundary Γ
3, it satisfies the Dirichlet boundary or Neumann boundary condition: (13)
u
1
0
,
t
=
θ
1
t
on
⁡
Γ
3
.On the boundary Γ
4, it satisfies the Dirichlet boundary conditions:(14)
u
2
d
,
t
=
θ
t
on
⁡
Γ
4
.On the phase change moving interface (Γ
5), it satisfies the condition of temperature continuity and the Stefan condition:(15)
u
1
ξ
t
,
t
=
u
2
ξ
t
,
t
=
u
∗
,where u
∗ is the phase change temperature. Hence,(16)
κ
d
ξ
t
t
=
k
2
∂
u
2
x
,
t
∂
x
x
=
ξ
t
-
k
1
∂
u
1
x
,
t
∂
x
x
=
ξ
t
,where κ is the latent heat of fusion per unit volume and k
1 and k
2 are the thermal conductivity in the solid and liquid phase, respectively. The final time internal temperature on the boundary Γ
6 satisfies(17)
u
1
x
,
t
∗
=
u
1
T
x
on
⁡
Γ
6
.


## 3. The Solution of Two-Phase Inverse Stefan Problem

### 3.1. The Homotopy Perturbation Method

The homotopy perturbation method was presented in [13, 14]; it arose as a combination of two other methods: the perturbation method and the homotopy technique from topology. Homotopy perturbation method was applied on differential equations, partial differential equations, and the direct and inverse problems of heat transfer by many authors [15–17].Homotopy perturbation method was used to find the exact or the approximate solutions of the linear and nonlinear integral equations [16, 18], the two-dimensional integral equations [19–21], and the integrodifferential equations [22–25]. This method enables seeking a solution of the following operator equation:(18)
B
w
-
ζ
z
=
0
z
∈
Ω
,where B is an operator, w is a sought function, and ζ is a known function on the domain of Ω. The operator B is presented in the form of sum:(19)
B
v
=
R
w
+
Q
w
,where R defines the linear operator and Q is the remaining part of operator B. So (18) can be written as(20)
R
w
+
Q
w
-
ζ
z
=
0
z
∈
Ω
.We define a new operator H which is called the homotopy operator in the following:(21)
H
v
,
p
=
1
-
p
R
v
-
R
w
0
+
p
B
v
-
ζ
z
=
0
,where p
∈
[
0,1
] is an embedding parameter, v
(
z
,
p
)
:
Ω
×
[
0,1
]
→
R, and w
0 defines the initial approximation of a solution of (18); (21) is transformed into the following form:(22)
H
v
,
p
=
R
v
-
R
w
0
+
p
R
w
0
+
p
Q
v
-
ζ
z
=
0
.Therefore,(23)
H
v
,
0
=
R
v
-
R
w
0
,
H
v
,
1
=
R
w
0
+
Q
v
-
ζ
z
.For p
=
0 the solution of operator equation H
(
v
,
0
)
=
0 is equivalent to solution of problem R
(
v
)
-
R
(
w
0
)
=
0, whereas, for p
=
1, the solution of operator equation H
(
v
,
1
)
=
0 is equivalent to solution of problem (20). Thus, a monotonous change of parameter p from zero to one is just that continuous change of the trivial problem R
(
v
)
-
R
(
w
0
)
=
0 to the original problem. In topology, this is called homotopy. According to the HPM, we assume that the solution of (21) and (22) can be written as a power series in p:(24)
v
=
∑
i
=
0
∞
p
i
w
i
=
w
0
+
p
w
1
+
p
2
w
2
+
⋯
;setting p
=
1, the approximate solution of (20) is obtained:(25)
w
=
lim
p
→
1
⁡
v
=
∑
i
=
0
∞
w
i
=
w
0
+
w
1
+
w
2
+
⋯
;in most cases, the series in (25) is convergent, and in [15] an additional remark on the convergence of this series can be found.


### 3.2. The Improved Adomian Decomposition Method

The Adomian decomposition method was presented by Adomian [26]. This method was used to solve a wide variety of linear and nonlinear problems. Several other researchers had developed a modification to the ADM [27]. The modification usually involves a slight change and it is aimed at improving the accuracy of the series solution.The Adomian decomposition method is able to solve a wide class of nonlinear operator equations [28] in the form(26)
F
u
=
η
,where F
:
H
→
G is a nonlinear operator, η is a known element from Hilbert space G and u is the sought element from Hilbert space H. The operator F
(
u
) can be written as (27)
F
u
=
L
u
+
R
u
+
N
u
,where L is an invertible linear operator, R is a linear operator, and N is a nonlinear operator. The solution of (26) admits the decomposition into an infinite series of components(28)
u
=
∑
i
=
0
∞
g
i
;the nonlinear term N is equated to an infinite series(29)
N
u
=
∑
i
=
0
∞
A
i
;
A
n is the Adomian polynomials, which can be determined by(30)
A
0
=
N
g
0
,
(31)
A
n
=
1
n
!
d
n
d
λ
n
N
∑
i
=
0
n
λ
i
g
i
λ
=
0
n
≥
1
.Substituting (28), (29), (30), and (31) into operator equation (26) and using the inverse operator L
-
1, the following equation can be obtained:(32)
g
0
=
g
∗
+
L
-
1
η
,
g
n
=
-
L
-
1
R
g
n
-
1
-
L
-
1
A
n
-
1
,
n
≥
1
,where g
∗ is the function dependent on the initial and boundary conditions.The improved Adomian decomposition method [27]: it is important to note that the improved Adomian decomposition method is based on the assumption that the function y
=
g
∗
+
L
-
1
(
η
) can be divided into two parts, namely, y
1 and y
2. Under this assumption, the following equation can be obtained:(33)
y
=
y
1
+
y
2
.Accordingly, a slight variation is proposed on the components g
0 and g
1. The suggestion is that only the part y
1 is assigned to the component g
0, whereas the remaining part y
2 is combined with other terms given in (33) to define g
1. Consequently, the recursive relation is written as follows:(34)
g
0
=
y
1
,
g
1
=
y
2
-
L
-
1
R
g
0
-
L
-
1
A
0
,
g
n
+
2
=
-
L
-
1
R
g
n
+
1
-
L
-
1
A
n
+
1
n
≥
0
.


### 3.3. The Solution of the Two-Phase Inverse Stefan Problem

In this section, we introduce the solution of the two-phase inverse Stefan problem. We split the two-phase inverse Stefan problem into two problems. The first is the determination ofu
1
(
x
,
t
) in domain D
1 and the moving interface ξ
(
t
) which is satisfying (7), (8), (10), and (13). Once the boundary x
=
ξ
(
t
) and the heat flux (
∂
u
2
(
x
,
t
)
/
∂
x
)
x
=
ξ
(
t
) have been obtained, the second problem for determining the temperature u
2
(
x
,
t
) in domain D
2 is solved.The solution ofu
1
(
x
,
t
) and moving interface ξ
(
t
) based on HPM: because the heat flux on the moving position is unknown in this two-phase inverse Stefan problem, the HPM is used to compute u
1
(
x
,
t
). In this section, the homotopy perturbation method is applied to solve the inverse problem. We make the homotopy map for (22):(35)
H
1
v
1
,
p
=
∂
2
v
1
∂
x
2
-
∂
2
u
1,0
∂
x
2
+
p
∂
2
u
1,0
∂
x
2
+
1
α
1
∂
v
i
∂
t
;the solution to this equation,(36)
H
1
v
1
,
p
=
0
,is sought in the form of a power series in p and it satisfies(37)
v
1
=
∑
i
=
0
∞
p
i
u
1
,
i
,where p
∈
[
0,1
]; by substituting (37) into (36), the following equation can be obtained:(38)
∑
i
=
0
∞
p
i
∂
2
u
1
,
i
∂
x
2
=
∂
2
u
1,0
∂
x
2
-
p
∂
2
u
1,0
∂
x
2
+
1
α
1
∑
i
=
1
∞
∂
u
1
,
i
-
1
∂
t
.Next, by comparing the same powers of parameter p in (38), the following equation can be obtained:(39)
∂
2
u
1,1
∂
x
2
=
1
α
1
∂
u
1,0
∂
t
-
∂
2
u
1,0
∂
x
2
,
(40)
∂
2
u
1
,
i
∂
x
2
=
1
α
1
∂
u
1
,
i
-
1
∂
t
for
i
≥
2
.Equations (39) and (40) must be supplemented by the boundary conditions. Therefore, for (39), we set the following boundary conditions:(41)
u
1,0
0
,
t
+
u
1,1
0
,
t
=
θ
1
t
,
(42)
u
1,0
ξ
t
,
t
+
u
1,1
ξ
t
,
t
=
u
∗
,and for (40) we set the conditions in the following (i
≥
2):(43)
u
1
,
i
0
,
t
=
0
,
(44)
u
1
,
i
ξ
t
,
t
=
0
.The above conditions are selected. So we give the n-order approximate solution: (45)
u
^
1
,
n
=
∑
i
=
1
n
u
1
,
i
,because u
1
,
i are determined by the unknown function ξ
(
t
). This function is derived in the form(46)
ξ
t
=
∑
i
=
1
m
c
i
ψ
i
t
.Thus we define a cost function J
1 according to (10) and (17):(47)
J
1
c
1
,
…
,
c
m
=
∫
0
ξ
0
u
1
,
n
x
,
0
-
φ
1
x
2
d
x
+
∫
0
ξ
t
∗
u
1
,
n
x
,
t
∗
-
u
1
T
x
2
d
x
.Therefore, c
i are chosen according to the minimum of (47). The moving interface ξ
(
t
) and the solution of u
1 in domain D
1 can be obtained.The boundary condition (13) on the boundary Γ
3 and the condition of temperature continuity (15) should be fulfilled for n
≥
1. In this way, while looking for the solution of this problem in domain D
1, the initial approximation u
1,1 will be found by solving the system of (39) with conditions (41) and (42). Conversely, functions u
1
,
i, i
=
1,2
,
…, are determined by the recurrent solution of (40) with conditions (42) and (44). Consequently, the moving interface ξ
(
t
) and the solution of u
1 in domain D
1 can be obtained by the homotopy perturbation method with optimization. According to the Stefan condition (16) on the phase change moving interface Γ
5, the temperature u
2
(
x
,
t
) in domain D
2 can be recovered by IADM with optimization method.The solution ofu
2
(
x
,
t
) based on IADM with optimization: in this problem, we consider the operator equations (27), where(48)
L
u
2
=
∂
2
u
2
∂
x
2
,
R
u
2
=
-
1
α
2
∂
u
2
∂
t
,
N
u
2
=
0
,
η
=
0
.The inverse operator L
-
1 is given by the following equation:(49)
L
-
1
u
2
=
∫
d
x
∫
ξ
t
x
u
2
x
,
t
d
x
d
x
.The boundary condition and the Stefan condition from (14) are used to obtain the following equation:(50)
L
-
1
L
∂
2
u
2
∂
x
2
=
u
2
x
,
t
-
θ
t
-
∂
u
2
x
,
t
∂
x
x
=
ξ
t
x
-
d
.Hence(51)
g
∗
=
θ
t
+
∂
u
2
x
,
t
∂
x
x
=
ξ
t
x
-
d
;then the following recursive relation can be obtained:(52)
g
0
=
u
2
d
,
t
=
θ
t
,
g
1
=
∂
u
2
x
,
t
∂
x
x
=
ξ
t
x
-
d
+
1
α
2
∫
d
x
∫
ξ
t
x
∂
g
0
∂
t
d
x
d
x
,
g
n
+
2
=
1
α
2
∫
d
x
∫
ξ
t
x
∂
g
n
-
1
∂
t
d
x
d
x
,
n
≥
0
;because u
1
(
x
,
t
) and ξ
(
t
) are determined by the above HPM with optimization method, (
∂
u
2
(
x
,
t
)
/
∂
x
)
x
=
ξ
(
t
) can be obtained by (16). An approximation solution is expressed in the form(53)
u
2
,
n
=
∑
i
=
0
n
g
i
n
∈
N
,and in this inverse Stefan problem, θ
(
t
) is unknown on the boundary Γ
4 and it is derived in the following form:(54)
θ
t
=
∑
i
=
1
m
q
i
ϕ
i
t
,where q
i
∈
R and ϕ
i
(
t
) are the basis function. The coefficients q
i are selected to show minimal deviation of function (55) from conditions (11) and (15). Thus q
i are chosen according to the minimum of the following functional:(55)
J
2
q
1
,
…
,
q
m
=
∫
0
t
∗
u
2
,
n
ξ
t
,
t
-
u
∗
2
d
t
+
∫
ξ
0
d
u
2
,
n
x
,
0
-
φ
2
x
2
d
x
,where J
2
(
q
1
,
…
,
q
m
) is a nonlinear function. The optimization method plays a very important role in solving this two-phase inverse Stefan problem. Therefore, a modified spectral DY conjugate gradient method is presented. In Section 4, this optimization method is introduced.Remark. The IADM with optimization method can be used to solve the direct and inverse one-phase Stefan problem which is designed in the literature [8, 9]. And this method is better than ADM with optimization method. In Section 5, we give two examples of one-phase Stefan problem to illustrate the effectiveness of the IADM with optimization method.


## 3.1. The Homotopy Perturbation Method

The homotopy perturbation method was presented in [13, 14]; it arose as a combination of two other methods: the perturbation method and the homotopy technique from topology. Homotopy perturbation method was applied on differential equations, partial differential equations, and the direct and inverse problems of heat transfer by many authors [15–17].Homotopy perturbation method was used to find the exact or the approximate solutions of the linear and nonlinear integral equations [16, 18], the two-dimensional integral equations [19–21], and the integrodifferential equations [22–25]. This method enables seeking a solution of the following operator equation:(18)
B
w
-
ζ
z
=
0
z
∈
Ω
,where B is an operator, w is a sought function, and ζ is a known function on the domain of Ω. The operator B is presented in the form of sum:(19)
B
v
=
R
w
+
Q
w
,where R defines the linear operator and Q is the remaining part of operator B. So (18) can be written as(20)
R
w
+
Q
w
-
ζ
z
=
0
z
∈
Ω
.We define a new operator H which is called the homotopy operator in the following:(21)
H
v
,
p
=
1
-
p
R
v
-
R
w
0
+
p
B
v
-
ζ
z
=
0
,where p
∈
[
0,1
] is an embedding parameter, v
(
z
,
p
)
:
Ω
×
[
0,1
]
→
R, and w
0 defines the initial approximation of a solution of (18); (21) is transformed into the following form:(22)
H
v
,
p
=
R
v
-
R
w
0
+
p
R
w
0
+
p
Q
v
-
ζ
z
=
0
.Therefore,(23)
H
v
,
0
=
R
v
-
R
w
0
,
H
v
,
1
=
R
w
0
+
Q
v
-
ζ
z
.For p
=
0 the solution of operator equation H
(
v
,
0
)
=
0 is equivalent to solution of problem R
(
v
)
-
R
(
w
0
)
=
0, whereas, for p
=
1, the solution of operator equation H
(
v
,
1
)
=
0 is equivalent to solution of problem (20). Thus, a monotonous change of parameter p from zero to one is just that continuous change of the trivial problem R
(
v
)
-
R
(
w
0
)
=
0 to the original problem. In topology, this is called homotopy. According to the HPM, we assume that the solution of (21) and (22) can be written as a power series in p:(24)
v
=
∑
i
=
0
∞
p
i
w
i
=
w
0
+
p
w
1
+
p
2
w
2
+
⋯
;setting p
=
1, the approximate solution of (20) is obtained:(25)
w
=
lim
p
→
1
⁡
v
=
∑
i
=
0
∞
w
i
=
w
0
+
w
1
+
w
2
+
⋯
;in most cases, the series in (25) is convergent, and in [15] an additional remark on the convergence of this series can be found.


## 3.2. The Improved Adomian Decomposition Method

The Adomian decomposition method was presented by Adomian [26]. This method was used to solve a wide variety of linear and nonlinear problems. Several other researchers had developed a modification to the ADM [27]. The modification usually involves a slight change and it is aimed at improving the accuracy of the series solution.The Adomian decomposition method is able to solve a wide class of nonlinear operator equations [28] in the form(26)
F
u
=
η
,where F
:
H
→
G is a nonlinear operator, η is a known element from Hilbert space G and u is the sought element from Hilbert space H. The operator F
(
u
) can be written as (27)
F
u
=
L
u
+
R
u
+
N
u
,where L is an invertible linear operator, R is a linear operator, and N is a nonlinear operator. The solution of (26) admits the decomposition into an infinite series of components(28)
u
=
∑
i
=
0
∞
g
i
;the nonlinear term N is equated to an infinite series(29)
N
u
=
∑
i
=
0
∞
A
i
;
A
n is the Adomian polynomials, which can be determined by(30)
A
0
=
N
g
0
,
(31)
A
n
=
1
n
!
d
n
d
λ
n
N
∑
i
=
0
n
λ
i
g
i
λ
=
0
n
≥
1
.Substituting (28), (29), (30), and (31) into operator equation (26) and using the inverse operator L
-
1, the following equation can be obtained:(32)
g
0
=
g
∗
+
L
-
1
η
,
g
n
=
-
L
-
1
R
g
n
-
1
-
L
-
1
A
n
-
1
,
n
≥
1
,where g
∗ is the function dependent on the initial and boundary conditions.The improved Adomian decomposition method [27]: it is important to note that the improved Adomian decomposition method is based on the assumption that the function y
=
g
∗
+
L
-
1
(
η
) can be divided into two parts, namely, y
1 and y
2. Under this assumption, the following equation can be obtained:(33)
y
=
y
1
+
y
2
.Accordingly, a slight variation is proposed on the components g
0 and g
1. The suggestion is that only the part y
1 is assigned to the component g
0, whereas the remaining part y
2 is combined with other terms given in (33) to define g
1. Consequently, the recursive relation is written as follows:(34)
g
0
=
y
1
,
g
1
=
y
2
-
L
-
1
R
g
0
-
L
-
1
A
0
,
g
n
+
2
=
-
L
-
1
R
g
n
+
1
-
L
-
1
A
n
+
1
n
≥
0
.


## 3.3. The Solution of the Two-Phase Inverse Stefan Problem

In this section, we introduce the solution of the two-phase inverse Stefan problem. We split the two-phase inverse Stefan problem into two problems. The first is the determination ofu
1
(
x
,
t
) in domain D
1 and the moving interface ξ
(
t
) which is satisfying (7), (8), (10), and (13). Once the boundary x
=
ξ
(
t
) and the heat flux (
∂
u
2
(
x
,
t
)
/
∂
x
)
x
=
ξ
(
t
) have been obtained, the second problem for determining the temperature u
2
(
x
,
t
) in domain D
2 is solved.The solution ofu
1
(
x
,
t
) and moving interface ξ
(
t
) based on HPM: because the heat flux on the moving position is unknown in this two-phase inverse Stefan problem, the HPM is used to compute u
1
(
x
,
t
). In this section, the homotopy perturbation method is applied to solve the inverse problem. We make the homotopy map for (22):(35)
H
1
v
1
,
p
=
∂
2
v
1
∂
x
2
-
∂
2
u
1,0
∂
x
2
+
p
∂
2
u
1,0
∂
x
2
+
1
α
1
∂
v
i
∂
t
;the solution to this equation,(36)
H
1
v
1
,
p
=
0
,is sought in the form of a power series in p and it satisfies(37)
v
1
=
∑
i
=
0
∞
p
i
u
1
,
i
,where p
∈
[
0,1
]; by substituting (37) into (36), the following equation can be obtained:(38)
∑
i
=
0
∞
p
i
∂
2
u
1
,
i
∂
x
2
=
∂
2
u
1,0
∂
x
2
-
p
∂
2
u
1,0
∂
x
2
+
1
α
1
∑
i
=
1
∞
∂
u
1
,
i
-
1
∂
t
.Next, by comparing the same powers of parameter p in (38), the following equation can be obtained:(39)
∂
2
u
1,1
∂
x
2
=
1
α
1
∂
u
1,0
∂
t
-
∂
2
u
1,0
∂
x
2
,
(40)
∂
2
u
1
,
i
∂
x
2
=
1
α
1
∂
u
1
,
i
-
1
∂
t
for
i
≥
2
.Equations (39) and (40) must be supplemented by the boundary conditions. Therefore, for (39), we set the following boundary conditions:(41)
u
1,0
0
,
t
+
u
1,1
0
,
t
=
θ
1
t
,
(42)
u
1,0
ξ
t
,
t
+
u
1,1
ξ
t
,
t
=
u
∗
,and for (40) we set the conditions in the following (i
≥
2):(43)
u
1
,
i
0
,
t
=
0
,
(44)
u
1
,
i
ξ
t
,
t
=
0
.The above conditions are selected. So we give the n-order approximate solution: (45)
u
^
1
,
n
=
∑
i
=
1
n
u
1
,
i
,because u
1
,
i are determined by the unknown function ξ
(
t
). This function is derived in the form(46)
ξ
t
=
∑
i
=
1
m
c
i
ψ
i
t
.Thus we define a cost function J
1 according to (10) and (17):(47)
J
1
c
1
,
…
,
c
m
=
∫
0
ξ
0
u
1
,
n
x
,
0
-
φ
1
x
2
d
x
+
∫
0
ξ
t
∗
u
1
,
n
x
,
t
∗
-
u
1
T
x
2
d
x
.Therefore, c
i are chosen according to the minimum of (47). The moving interface ξ
(
t
) and the solution of u
1 in domain D
1 can be obtained.The boundary condition (13) on the boundary Γ
3 and the condition of temperature continuity (15) should be fulfilled for n
≥
1. In this way, while looking for the solution of this problem in domain D
1, the initial approximation u
1,1 will be found by solving the system of (39) with conditions (41) and (42). Conversely, functions u
1
,
i, i
=
1,2
,
…, are determined by the recurrent solution of (40) with conditions (42) and (44). Consequently, the moving interface ξ
(
t
) and the solution of u
1 in domain D
1 can be obtained by the homotopy perturbation method with optimization. According to the Stefan condition (16) on the phase change moving interface Γ
5, the temperature u
2
(
x
,
t
) in domain D
2 can be recovered by IADM with optimization method.The solution ofu
2
(
x
,
t
) based on IADM with optimization: in this problem, we consider the operator equations (27), where(48)
L
u
2
=
∂
2
u
2
∂
x
2
,
R
u
2
=
-
1
α
2
∂
u
2
∂
t
,
N
u
2
=
0
,
η
=
0
.The inverse operator L
-
1 is given by the following equation:(49)
L
-
1
u
2
=
∫
d
x
∫
ξ
t
x
u
2
x
,
t
d
x
d
x
.The boundary condition and the Stefan condition from (14) are used to obtain the following equation:(50)
L
-
1
L
∂
2
u
2
∂
x
2
=
u
2
x
,
t
-
θ
t
-
∂
u
2
x
,
t
∂
x
x
=
ξ
t
x
-
d
.Hence(51)
g
∗
=
θ
t
+
∂
u
2
x
,
t
∂
x
x
=
ξ
t
x
-
d
;then the following recursive relation can be obtained:(52)
g
0
=
u
2
d
,
t
=
θ
t
,
g
1
=
∂
u
2
x
,
t
∂
x
x
=
ξ
t
x
-
d
+
1
α
2
∫
d
x
∫
ξ
t
x
∂
g
0
∂
t
d
x
d
x
,
g
n
+
2
=
1
α
2
∫
d
x
∫
ξ
t
x
∂
g
n
-
1
∂
t
d
x
d
x
,
n
≥
0
;because u
1
(
x
,
t
) and ξ
(
t
) are determined by the above HPM with optimization method, (
∂
u
2
(
x
,
t
)
/
∂
x
)
x
=
ξ
(
t
) can be obtained by (16). An approximation solution is expressed in the form(53)
u
2
,
n
=
∑
i
=
0
n
g
i
n
∈
N
,and in this inverse Stefan problem, θ
(
t
) is unknown on the boundary Γ
4 and it is derived in the following form:(54)
θ
t
=
∑
i
=
1
m
q
i
ϕ
i
t
,where q
i
∈
R and ϕ
i
(
t
) are the basis function. The coefficients q
i are selected to show minimal deviation of function (55) from conditions (11) and (15). Thus q
i are chosen according to the minimum of the following functional:(55)
J
2
q
1
,
…
,
q
m
=
∫
0
t
∗
u
2
,
n
ξ
t
,
t
-
u
∗
2
d
t
+
∫
ξ
0
d
u
2
,
n
x
,
0
-
φ
2
x
2
d
x
,where J
2
(
q
1
,
…
,
q
m
) is a nonlinear function. The optimization method plays a very important role in solving this two-phase inverse Stefan problem. Therefore, a modified spectral DY conjugate gradient method is presented. In Section 4, this optimization method is introduced.Remark. The IADM with optimization method can be used to solve the direct and inverse one-phase Stefan problem which is designed in the literature [8, 9]. And this method is better than ADM with optimization method. In Section 5, we give two examples of one-phase Stefan problem to illustrate the effectiveness of the IADM with optimization method.


## 4. A Modified Spectral DY Conjugate Gradient Method

Conjugate gradient methods are very efficient for solving the unconstrained optimization problemm
i
n
x
∈
R
n
f
(
x
). Their iterative formula is given as follows:(56)
x
k
+
1
=
x
k
+
α
k
d
k
,
d
k
=
-
g
k
,
if
k
=
1
;
-
g
k
+
β
k
d
k
-
1
if
k
≥
1
,where step-size α
k is positive, g
k
=
∇
f
(
x
k
), and β
k is a scalar. In addition, α
k is a step length which is calculated by the line search. β
k is a very important factor in conjugate gradient methods; the FR, PRP, DY, and CD methods [29] are used for choosing this scalar. Among them, the DY method is regarded as one of the efficient conjugate gradient methods. The scalar β
k is given by(57)
β
k
D
Y
=
g
k
d
k
-
1
T
y
k
-
1
,where y
k
-
1
=
g
k
-
g
k
-
1 and · stands for Euclidean norm. Recently, Birgin and Martłnez [30] presented a spectral conjugate gradient method which is combining conjugate gradient method and spectral gradient method. And the direction d
k is obtained by(58)
d
k
=
-
θ
k
g
k
+
β
k
d
k
-
1
,
β
k
=
θ
k
y
k
-
1
-
s
k
-
1
T
g
k
d
k
-
1
T
y
k
-
1
,where θ
k
=
s
k
-
1
T
s
k
-
1
/
s
k
-
1
T
y
k
-
1 is a parameter with s
k
-
1
=
x
k
-
x
k
-
1 and θ
k is regarded as spectral gradient. However, according to [31], the spectral conjugate gradient method can not guarantee producing the descent directions.Based on the above analysis, in the following, we describe a modified spectral DY conjugate gradient method and give the convergence of this method with the Wolfe type line search rules.


### 4.1. The Modified Spectral DY Conjugate Gradient Method

The new method is introduced:(59)
d
k
=
-
g
k
,
if
k
=
1
;
-
θ
k
g
k
+
β
k
d
k
-
1
if
k
≥
1
,where θ
k
=
1
+
g
k
-
1
2
/
g
k
2 and β
k
=
β
k
D
Y. This method is called the MDY method.Algorithm 1.
Consider the following.Step 1.
The following are given:x
0
∈
R
n, ε
>
0, and the max iteration step N
m
a
x.Step 2.
Setk
=
1; compute d
1
=
-
g
1
=
-
∇
f
(
x
0
); if g
1
≤
ε, then stop.Step 3.
Findα
k according to Wolfe type line search rule.Step 4.
Letx
k
+
1
=
x
k
+
α
k
d
k and g
k
+
1
=
g
(
x
k
+
1
); if g
k
+
1
≤
ε, then stop.Step 5.
Computed
k
+
1 by (59).Step 6.
Setk
=
k
+
1; go to Step 3.


### 4.2. Convergence Analysis

From the above, the MDY conjugate gradient method provides a descent direction for the objective function. In the following, we give the convergence of this algorithm with the Wolfe type line search rules. Firstly, we introduce some assumptions:(i)
The level setΩ
=
{
x
∈
R
n
∣
f
(
x
)
≤
f
(
x
1
)
} is bounded.
(ii)
In some neighborhoodN of Ω, f
(
x
) is continuously differentiable and its gradient is Lipschitz continuous; namely, there exists a constant L
>
0 such that(60)
g
x
-
g
y
≤
L
x
-
y
,
∀
x
,
y
∈
N
.It is implied from the above assumptions that there exist two positive constants β and γ such that(61)
x
≤
B
,
g
x
≤
γ
,
∀
x
∈
Ω
.The Wolfe type line search [29] is(62)
f
x
k
-
f
x
k
+
α
k
d
k
≥
δ
α
k
2
d
k
2
,
(63)
g
x
k
+
α
k
d
k
T
d
k
≥
-
2
δ
α
k
d
k
2
,and the strong Wolfe line search is(64)
f
x
k
+
α
k
d
k
≤
f
x
k
+
δ
α
k
g
k
T
d
k
,
(65)
g
x
k
+
α
k
d
k
T
d
k
≤
σ
g
k
T
d
k
,where 0
<
δ
<
σ
<
1.Lemma 2.
Let the sequences{
g
k
} and {
d
k
} be generated by the proposed algorithm (Algorithm 1) with Wolfe type line search rules; then g
k
T
d
k
<
0 and(66)
β
k
D
Y
≤
g
k
T
d
k
g
k
-
1
T
d
k
-
1
.Proof.
According to (59), we have(67)
g
k
T
d
k
=
-
θ
k
g
k
2
+
β
k
D
Y
g
k
T
d
k
-
1
=
β
k
D
Y
-
θ
k
d
k
-
1
g
k
T
-
g
k
-
1
T
+
g
k
T
d
k
-
1
=
β
k
D
Y
g
k
-
1
T
d
k
-
1
-
g
k
-
1
2
d
k
-
1
y
k
-
1
T
g
k
2
≤
β
k
D
Y
g
k
-
1
T
d
k
-
1
.If we choose k
=
1, then(68)
g
1
T
d
1
=
-
g
1
2
<
0
;according to (65), we have(69)
y
1
T
d
1
≥
σ
-
1
g
1
T
d
1
≥
0
.So β
2
D
Y
≥
0. For k
>
1, according to (67), we have(70)
g
k
T
d
k
<
0
,by induction. So we obtain(71)
β
k
D
Y
≤
g
k
T
d
k
g
k
-
1
T
d
k
-
1
;the proof is completed.Theorem 3.
Let the sequences{
g
k
} and {
d
k
} be generated by the proposed algorithm (Algorithm 1) with Wolfe type line search rules; then(72)
g
k
T
d
k
≤
-
C
¯
∥
g
k
∥
2
,where C
¯
>
0 is a constant.Proof.
According to (59), we have(73)
g
k
T
d
k
=
-
θ
k
g
k
2
+
β
k
D
Y
g
k
T
d
k
-
1
=
-
g
k
-
1
2
+
g
k
T
d
k
-
1
d
k
-
1
T
y
k
-
1
-
1
g
k
2
≤
g
k
-
1
T
d
k
-
1
d
k
-
1
T
y
k
-
1
g
k
2
.From the Wolfe type line search rules (65), we can obtain(74)
-
σ
+
1
g
k
-
1
T
d
k
-
1
≥
d
k
-
1
T
g
k
-
g
k
-
1
>
0
;according to (73) and (74), we have (75)
g
k
T
d
k
≤
-
1
σ
+
1
g
k
2
;we define C
¯
=
1
/
(
σ
+
1
). Therefore, we can obtain inequality (72).Lemma 4.
Suppose that assumptions (i) and (ii) hold; let{
g
k
} and {
d
k
} be generated by the proposed algorithm (Algorithm 1) with Wolfe type line search rules (62) and (63); then one has(76)
∑
k
≥
0
g
k
T
d
k
2
d
k
2
≤
∞
.Proof.
From assumption (ii), we have(77)
g
x
k
+
α
k
d
k
T
-
g
x
k
T
≤
L
α
k
d
k
,so we can obtain(78)
g
x
k
+
α
k
d
k
T
d
k
-
g
x
k
T
d
k
≤
L
α
k
d
k
2
.From inequality (63), we obtain(79)
2
σ
α
k
d
k
2
≥
-
g
x
k
+
α
k
d
k
T
d
k
.According to inequalities (78) and (79), we have(80)
2
σ
+
L
α
k
d
k
2
≥
-
g
k
T
d
k
.So we can obtain(81)
α
k
d
k
≥
1
2
σ
+
L
-
g
k
T
d
k
d
k
.Owing to inequality (81), we have(82)
∑
k
=
1
∞
-
g
k
T
d
k
d
k
≤
2
σ
+
L
2
∑
k
=
1
∞
α
k
2
d
k
2
≤
2
σ
+
L
2
∑
k
=
1
∞
f
x
k
-
f
x
k
+
1
<
+
∞
.From (82), we can obtain inequality (76).Theorem 5.
Suppose that assumptions (i) and (ii) hold. Let{
g
k
} and {
d
k
} be generated by the algorithm with Wolfe type line search rules. Then one has(83)
liminf
k
→
∞
⁡
g
k
=
0
.Proof.
By contradiction, we suppose that there exists a positive constantγ
¯
>
0, such that(84)
g
k
≥
γ
¯holds for all k
>
1. From (59), we have(85)
d
k
2
=
-
2
θ
k
d
k
T
g
k
+
β
k
D
Y
2
d
k
-
1
2
-
θ
k
2
g
k
2
.According to Lemma 4 and the above equation, we have(86)
d
k
2
≤
g
k
T
d
k
g
k
-
1
T
d
k
-
1
2
d
k
-
1
2
-
θ
k
2
g
k
2
-
2
θ
k
d
k
T
g
k
.Multiplying the above inequality by 1
/
g
k
T
d
k
2, we have(87)
d
k
2
g
k
T
d
k
2
≤
d
k
-
1
2
g
k
-
1
T
d
k
-
1
2
-
θ
k
g
k
g
k
T
d
k
+
1
g
k
2
+
1
g
k
2
≤
d
k
-
1
2
g
k
-
1
T
d
k
-
1
2
+
1
g
k
2
≤
∑
i
=
1
k
1
g
i
2
2
≤
1
γ
¯
2
k
.Then from inequality (83), we have(88)
∑
k
=
1
∞
g
k
T
d
k
2
g
k
2
≥
γ
¯
2
∑
k
=
1
∞
1
k
=
∞
,which contradicts (76). Therefore, (83) holds.


## 4.1. The Modified Spectral DY Conjugate Gradient Method

The new method is introduced:(59)
d
k
=
-
g
k
,
if
k
=
1
;
-
θ
k
g
k
+
β
k
d
k
-
1
if
k
≥
1
,where θ
k
=
1
+
g
k
-
1
2
/
g
k
2 and β
k
=
β
k
D
Y. This method is called the MDY method.Algorithm 1.
Consider the following.Step 1.
The following are given:x
0
∈
R
n, ε
>
0, and the max iteration step N
m
a
x.Step 2.
Setk
=
1; compute d
1
=
-
g
1
=
-
∇
f
(
x
0
); if g
1
≤
ε, then stop.Step 3.
Findα
k according to Wolfe type line search rule.Step 4.
Letx
k
+
1
=
x
k
+
α
k
d
k and g
k
+
1
=
g
(
x
k
+
1
); if g
k
+
1
≤
ε, then stop.Step 5.
Computed
k
+
1 by (59).Step 6.
Setk
=
k
+
1; go to Step 3.


## 4.2. Convergence Analysis

From the above, the MDY conjugate gradient method provides a descent direction for the objective function. In the following, we give the convergence of this algorithm with the Wolfe type line search rules. Firstly, we introduce some assumptions:(i)
The level setΩ
=
{
x
∈
R
n
∣
f
(
x
)
≤
f
(
x
1
)
} is bounded.
(ii)
In some neighborhoodN of Ω, f
(
x
) is continuously differentiable and its gradient is Lipschitz continuous; namely, there exists a constant L
>
0 such that(60)
g
x
-
g
y
≤
L
x
-
y
,
∀
x
,
y
∈
N
.It is implied from the above assumptions that there exist two positive constants β and γ such that(61)
x
≤
B
,
g
x
≤
γ
,
∀
x
∈
Ω
.The Wolfe type line search [29] is(62)
f
x
k
-
f
x
k
+
α
k
d
k
≥
δ
α
k
2
d
k
2
,
(63)
g
x
k
+
α
k
d
k
T
d
k
≥
-
2
δ
α
k
d
k
2
,and the strong Wolfe line search is(64)
f
x
k
+
α
k
d
k
≤
f
x
k
+
δ
α
k
g
k
T
d
k
,
(65)
g
x
k
+
α
k
d
k
T
d
k
≤
σ
g
k
T
d
k
,where 0
<
δ
<
σ
<
1.Lemma 2.
Let the sequences{
g
k
} and {
d
k
} be generated by the proposed algorithm (Algorithm 1) with Wolfe type line search rules; then g
k
T
d
k
<
0 and(66)
β
k
D
Y
≤
g
k
T
d
k
g
k
-
1
T
d
k
-
1
.Proof.
According to (59), we have(67)
g
k
T
d
k
=
-
θ
k
g
k
2
+
β
k
D
Y
g
k
T
d
k
-
1
=
β
k
D
Y
-
θ
k
d
k
-
1
g
k
T
-
g
k
-
1
T
+
g
k
T
d
k
-
1
=
β
k
D
Y
g
k
-
1
T
d
k
-
1
-
g
k
-
1
2
d
k
-
1
y
k
-
1
T
g
k
2
≤
β
k
D
Y
g
k
-
1
T
d
k
-
1
.If we choose k
=
1, then(68)
g
1
T
d
1
=
-
g
1
2
<
0
;according to (65), we have(69)
y
1
T
d
1
≥
σ
-
1
g
1
T
d
1
≥
0
.So β
2
D
Y
≥
0. For k
>
1, according to (67), we have(70)
g
k
T
d
k
<
0
,by induction. So we obtain(71)
β
k
D
Y
≤
g
k
T
d
k
g
k
-
1
T
d
k
-
1
;the proof is completed.Theorem 3.
Let the sequences{
g
k
} and {
d
k
} be generated by the proposed algorithm (Algorithm 1) with Wolfe type line search rules; then(72)
g
k
T
d
k
≤
-
C
¯
∥
g
k
∥
2
,where C
¯
>
0 is a constant.Proof.
According to (59), we have(73)
g
k
T
d
k
=
-
θ
k
g
k
2
+
β
k
D
Y
g
k
T
d
k
-
1
=
-
g
k
-
1
2
+
g
k
T
d
k
-
1
d
k
-
1
T
y
k
-
1
-
1
g
k
2
≤
g
k
-
1
T
d
k
-
1
d
k
-
1
T
y
k
-
1
g
k
2
.From the Wolfe type line search rules (65), we can obtain(74)
-
σ
+
1
g
k
-
1
T
d
k
-
1
≥
d
k
-
1
T
g
k
-
g
k
-
1
>
0
;according to (73) and (74), we have (75)
g
k
T
d
k
≤
-
1
σ
+
1
g
k
2
;we define C
¯
=
1
/
(
σ
+
1
). Therefore, we can obtain inequality (72).Lemma 4.
Suppose that assumptions (i) and (ii) hold; let{
g
k
} and {
d
k
} be generated by the proposed algorithm (Algorithm 1) with Wolfe type line search rules (62) and (63); then one has(76)
∑
k
≥
0
g
k
T
d
k
2
d
k
2
≤
∞
.Proof.
From assumption (ii), we have(77)
g
x
k
+
α
k
d
k
T
-
g
x
k
T
≤
L
α
k
d
k
,so we can obtain(78)
g
x
k
+
α
k
d
k
T
d
k
-
g
x
k
T
d
k
≤
L
α
k
d
k
2
.From inequality (63), we obtain(79)
2
σ
α
k
d
k
2
≥
-
g
x
k
+
α
k
d
k
T
d
k
.According to inequalities (78) and (79), we have(80)
2
σ
+
L
α
k
d
k
2
≥
-
g
k
T
d
k
.So we can obtain(81)
α
k
d
k
≥
1
2
σ
+
L
-
g
k
T
d
k
d
k
.Owing to inequality (81), we have(82)
∑
k
=
1
∞
-
g
k
T
d
k
d
k
≤
2
σ
+
L
2
∑
k
=
1
∞
α
k
2
d
k
2
≤
2
σ
+
L
2
∑
k
=
1
∞
f
x
k
-
f
x
k
+
1
<
+
∞
.From (82), we can obtain inequality (76).Theorem 5.
Suppose that assumptions (i) and (ii) hold. Let{
g
k
} and {
d
k
} be generated by the algorithm with Wolfe type line search rules. Then one has(83)
liminf
k
→
∞
⁡
g
k
=
0
.Proof.
By contradiction, we suppose that there exists a positive constantγ
¯
>
0, such that(84)
g
k
≥
γ
¯holds for all k
>
1. From (59), we have(85)
d
k
2
=
-
2
θ
k
d
k
T
g
k
+
β
k
D
Y
2
d
k
-
1
2
-
θ
k
2
g
k
2
.According to Lemma 4 and the above equation, we have(86)
d
k
2
≤
g
k
T
d
k
g
k
-
1
T
d
k
-
1
2
d
k
-
1
2
-
θ
k
2
g
k
2
-
2
θ
k
d
k
T
g
k
.Multiplying the above inequality by 1
/
g
k
T
d
k
2, we have(87)
d
k
2
g
k
T
d
k
2
≤
d
k
-
1
2
g
k
-
1
T
d
k
-
1
2
-
θ
k
g
k
g
k
T
d
k
+
1
g
k
2
+
1
g
k
2
≤
d
k
-
1
2
g
k
-
1
T
d
k
-
1
2
+
1
g
k
2
≤
∑
i
=
1
k
1
g
i
2
2
≤
1
γ
¯
2
k
.Then from inequality (83), we have(88)
∑
k
=
1
∞
g
k
T
d
k
2
g
k
2
≥
γ
¯
2
∑
k
=
1
∞
1
k
=
∞
,which contradicts (76). Therefore, (83) holds.


## 5. Numerical Simulation Experiments

In order to illustrate the validity of the IADM with optimization method, we introduce the one-phase inverse Stefan problem in Example1. From the literature [9], we know the inverse Stefan problem consists in the calculation of the temperature distribution in the domain and in the reconstruction of the temperature distribution on the boundary, when the position of the function which describes the moving interface is known. And then, in Example 2, we solve the two-phase inverse Stefan problem based on the hybrid method which combines the Adomian decomposition method with the homotopy perturbation method.Example 1.
In this example, the one-phase inverse Stefan problem [9] is taken for the following values of parameters: α
1
=
1, k
1
=
1, u
∗
=
1, t
∗
=
1
/
2, κ
=
1
/
α
1, φ
1
(
x
)
=
exp
⁡
(
α
1
x
), and u
1
T
(
x
)
=
exp
⁡
(
α
1
/
2
-
x
). Then the exact solution of such formulated one-phase Stefan problem can be found from the following functions:(89)
u
1
x
,
t
=
exp
⁡
α
1
t
-
x
,
θ
1
t
=
exp
⁡
α
1
t
,
ξ
t
=
α
1
t
.The following equation,(90)
ϕ
i
t
=
exp
⁡
-
i
-
1
t
i
=
1,2
,
…
,
m
,is taken as basis function. The approximate solutions are compared with the exact solution and the values of the absolute errors are calculated from(91)
δ
h
=
∫
0
t
∗
h
e
t
-
h
r
t
d
t
1
/
2
,
δ
u
=
∫
∫
D
u
e
x
,
t
-
u
n
x
,
t
d
x
d
t
1
/
2
,where h
(
t
)
∈
{
θ
(
t
)
,
ξ
(
t
)
} (θ
(
t
) stands for the boundary condition on Γ
1 in the literature [9] and ξ
(
t
) is the moving interface), D is the domain of the temperature distribution, h
e
(
t
) is the exact value of function h
(
t
), h
r
(
t
) is the approximate value of function h
(
t
), u
e
(
x
,
t
) is the exact distribution of the temperature in the domain, and u
n
(
x
,
t
) is the reconstructed distribution of temperature. Percentage relative errors are calculated from(92)
Δ
h
=
δ
h
∫
0
t
∗
h
e
t
2
d
t
×
100
%
,
Δ
u
=
δ
u
∫
∫
D
u
e
x
,
t
2
d
x
d
t
×
100
%
.The comparison of the ADM and the IADM for reconstructing distribution of the temperature on boundaryΓ
1 in the literature [9] is shown in Figure 3 for n
=
1, n
=
2 and m
=
2, m
=
3. From Figure 3, the method based on the improved Adomian decomposition method (IADM) is better than the Adomian decomposition method (ADM).Figure 3
The reconstructing distribution of the temperature on boundary (Γ
1 in the literature [9]) for n
=
1, n
=
2 and m
=
2, m
=
3.
(a)
(b)
(c)Example 2.
This example is introduced to illustrate the effectiveness of the hybrid method which combines the Adomian decomposition method with the homotopy perturbation method, which is presented in the previous sections. Here, we consider an example of the two-phase inverse Stefan problem [12], in which α
1
=
5
/
2, α
2
=
5
/
4, k
1
=
6, k
2
=
2, κ
=
4
/
5, u
∗
=
1, t
∗
=
1, d
=
3, s
=
3
/
2, and φ
1
(
x
)
=
e
(
3
-
2
x
)
/
10, φ
2
(
x
)
=
e
(
3
-
2
x
)
/
5, θ
1
(
t
)
=
e
(
t
+
3
)
/
10, and u
1
T
(
x
)
=
e
(
4
-
2
x
)
/
10. The exact solution of this inverse Stefan problem is described by the following functions:(93)
u
1
x
,
t
=
e
t
-
2
x
+
3
/
10
,
u
2
x
,
t
=
e
t
-
2
x
+
3
/
5
,
θ
t
=
e
t
-
3
/
5
,
ξ
t
=
t
+
3
2
.First, we use HPM with optimization method to solveu
1
(
x
,
t
) and ξ
(
t
). As the initial approximation u
1,0, this fulfills the initial condition:(94)
u
1,0
x
,
t
=
e
3
-
2
x
/
10
.By using the function u
1,0 in (39) and the conditions (41) and (42), the following system of equations can be obtained:(95)
∂
2
u
1,1
x
,
t
∂
x
2
=
-
1
25
e
3
-
2
x
/
10
.With the conditions(96)
u
1,1
0
,
t
=
e
3
/
10
e
t
/
10
-
1
,
u
1,1
x
,
t
x
=
ξ
t
=
1
-
e
3
-
2
ξ
t
/
10
,this system of equations can be solved. And the solution of this system is shown in the following:(97)
u
1,1
x
,
t
=
-
e
3
/
10
-
x
/
5
+
1
-
e
3
+
t
/
10
ξ
t
x
+
e
t
+
3
/
10
.Next, we designate the functionsu
1
,
i for i
≥
2 in the following:(98)
∂
2
u
1
,
i
∂
x
2
=
2
5
∂
u
1
,
i
-
1
∂
t
,and the conditions are shown:(99)
u
1
,
i
0
,
t
=
0
,
u
1
,
i
x
,
t
x
=
ξ
t
=
0
.By using (98) and the boundary conditions (99), we derive u
1
,
i
(
x
,
t
) for i
≥
2; according to (45), the approximate solution u
1
(
x
,
t
) is obtained. Here, we define basis functions according to (46):(100)
ψ
i
t
=
t
i
-
1
i
=
1
,
…
,
m
.We use the modified spectral DY (MDY) conjugate gradient method to solve function (47). The solution of moving interface ξ
(
t
) is shown in Figure 4. So the temperature u
1
(
x
,
t
) is derived in domain D
1. Values of the error in reconstruction of the position of the moving interface ξ
(
t
) and the temperature distribution u
1
(
x
,
t
) are shown in Table 1.Table 1
Values of the error in reconstruction of the position of the moving interface and the temperature distributionu
1
x
,
t.
HPM
δ
ξ
Δ
ξ
δ
u
1
δ
u
1
X
Δ
u
1
n
=
2
,
m
=
2
0.108712413045
0.61599738364
0.116978666115
0.109525848228
n
=
2
,
m
=
3
0.210642147214
1.19356205914
0.299360015653
0.280470260729
n
=
2
,
m
=
3
0.015061910569
0.08534533677
0.102526697222
0.096057215386Figure 4
The position of moving interfaceξ
(
t
).
(a)
(b)
(c)In order to verify the validity of the modified spectral DY (MDY) conjugate gradient method, we choosen
=
2, m
=
2, ε
=
1
0
-
4, the initial point c
0
=
(
2
.
2,2
.
5
), and the max step number N
m
a
x
=
100. The comparison results of DY and MDY are shown in Figure 5 and Table 2. Comparing to DY conjugate gradient method, the MDY conjugate gradient method is slightly better.Table 2
The comparison results of DY and MDY.
Algorithm
Iteration numbers
Optimization value
DY
50
(
1.52,0.48
)
MDY
23
(
1.51,0.501
)Figure 5
The cost functionJ
1 changes with the increasing of iteration numbers.The moving interfaceξ
(
t
) and u
1
(
x
,
t
) in domain D
1 are obtained by the homotopy perturbation method with optimization. Next, with the help of condition (16), we can determine the temperature u
2
(
x
,
t
) in domain D
2 by the improved Adomian decomposition method (IADM) with optimization. Without loss of generality, we choose the basis functions(101)
ϕ
i
t
=
t
i
-
1
,
i
=
1
,
…
,
m
.The exact and the reconstructed distribution of the temperature on boundaryΓ
4 are shown in Figure 6 for different number of basis functions ϕ
i
(
t
). And the errors in reconstruction of the temperature θ
(
t
) (on boundary Γ
4 in Figure 2) and the temperature distribution u
2
(
x
,
t
) are shown in Table 3. The inversion results of the function θ
(
t
) match the exact value very well.Table 3
Values of the error in reconstruction of the temperatureθ
(
t
) on boundary Γ
4 in Figure 2 and the temperature distribution u
2
x
,
t.
IADM
δ
θ
(
t
)
Δ
θ
(
t
)
δ
u
2
(
x
,
t
)
Δ
u
2
(
x
,
t
)
δ
u
2
X
(
x
,
t
)
n
=
2
,
m
=
2
0.0125536846
0.20525000280
0.649823007
0.5444948151
n
=
2
,
m
=
3
0.00907338187
0.14834781217
0.717820233
0.6014705398Figure 6
The reconstructing distribution of the temperatureθ
(
t
) on boundary Γ
4 in Figure 2.The measured temperatureu
1
T
(
x
) at the final time t
∗ on the boundary Γ
6 in Figure 2 is perturbed by the random error. So in the next section we consider the influence of measurement errors on the results. We set the measured temperature u
1
T
(
x
) with an error of 1
% and 3
%. We can obtain the moving interface ξ
(
t
) and θ
(
t
); the results are shown in Figures 7 and 8.Figure 7
The reconstructed position of moving interfaceξ
(
t
).Figure 8
The reconstructing distribution of the temperatureθ
(
t
) on boundary Γ
4 in Figure 2.In Figures7 and 8, the exact and the reconstructed distribution of the temperature on moving interface ξ
(
t
) and the boundary Γ
4 are shown, and the errors in reconstruction of the function describing the moving interface and boundary conditions are compiled in Table 4. The obtained reconstruction results illustrate that the functions ξ
(
t
) and θ
(
t
) are reconstructed very well. From the results, the hybrid method with optimization is stable with the input errors. When the input data are contaminated by the errors, the error of the reconstructed boundary conditions does not exceed the error of the input data.Table 4
Values of the error in reconstruction the moving interfaceξ
(
t
) and the temperature θ
(
t
) on boundary Γ
4 in Figure 2.
Perturbation error on theu
1
T
(
x
)
δ
ξ
(
t
)
Δ
ξ
(
t
)
δ
θ
(
t
)
Δ
θ
(
t
)
1% error
0.0968150593
0.771949532
0.032490326
0.7475026727
3% error
0.3725894970
2.9708217908
0.061914613
1.4244652005


## 6. Conclusion

Based on the homotopy perturbation method (HPM) with the improved Adomian decomposition method (IADM), a hybrid method with optimization is proposed to solve the two-phase inverse Stefan problem. This problem is more realistic from the practical point of view. The simulation experiment is used to verify this method. Experimental results match the exact value very well. Furthermore, this investigated method can also be used for solving one-phase direct Stefan problems.Optimization method is very crucial in this paper, so a modified spectral DY conjugate gradient method is presented. And we give the convergence of this MDY method. Simulation experiment illustrates the validity of this optimization method.



---

*Source: 319054-2015-06-09.xml*