One-dimensional heat equation with discontinuous conductance arXiv:1312.7396v1 [math.PR] 28 Dec 2013 Zhen-Qing Chen∗ and Mounir Zili August 5, 2021 Abstract We study a second-order parabolic equation with divergence form elliptic operator, having piecewise constant diffusion coefficients with two points of discontinuity. Such partial differential equations appear in the modelization of diffusion phenomena in medium consisting of three kind of materials. Using probabilistic methods, we present an explicit expression of the fundamental solution under certain conditions. We also derive small-time asymptotic expansion of the PDE’s solutions in the general case. The obtained results are directly usable in applications. AMS 2000 Mathematics Subject Classification: Primary 60H10, 60J55; Secondary 35K10, 35K08 Keywords and phrases: Stochastic differential equations, semimartingale local time, strong solution, heat kernel, asymptotic expansion 1 Introduction In many practical applications, one is encountered with heat propagations in heterogeneous media consisting of different kind of materials. Mathematically, such kind of heat propagation is modeled by heat equations with divergence form elliptic operators having discontinuous diffusion coefficients. In this paper we exam one-dimensional model where the medium consisting of three different kind of materials. Let a, ai , ρi , i = 1, 2, 3 be positive constants and define A(x) = a1 1{x≤0} + a2 1{0<x≤a} + a3 1{a<x} Set 1 d L= 2ρ(x) dx and  ρ(x) = ρ1 1{x≤0} + ρ2 1{0<x≤a} + ρ3 1{a<x} . d ρ(x)A(x) dx  , (1.1) which is a self-adjoint operator in L2 (R; ρ(x)dx). The parabolic equation ∂u(t, x) = Lu(t, x) ∂t ∗ Research partially supported by NSF Grant DMS-1206276 and NNSFC Grant 11128101. 1 (1.2) describes the heat propagation in the line that consists of three different kind of material. It is wellknown that there is an m-symmetric diffusion process X associated with L, where m(dx) := ρ(x)dx. Since both A(x) and ρ(x) are bounded between two positive constants, it is known (see [7]) that X possesses a jointly continuous transition density function q(t, x, y) with respect to the measure m(dx) on R that is symmetric in x and y and enjoys the following Aronson type Gaussian estimates: there exists positive constants C1 , C2 ≥ 1 so that for every t > 0 and x, y ∈ R, C1−1 t−1/2 exp(−C2 |x − y|2 /t) ≤ q(t, x, y) ≤ C1 t−1/2 exp(−|x − y|2 /(C2 t)). (1.3) Using Dirichlet form theory, the following semimartingale representation of X is derived in [10, Proposition 5]: p ρ2 a2 − ρ1 a1 b 0 ρ3 a3 − ρ2 a2 b a dXt = A(Xt )dBt + dLt + dLt , (1.4) ρ 2 a2 + ρ 1 a1 ρ3 a3 + ρ2 a2 b 0 and L b a are symmetric semimartingale local times where B is one-dimensional Brownian motion, L of X at 0 and a, respectively. Here for a semimartingale Y , its symmetric semimartingale local b w (Y ) at w ∈ R is defined to be time L t Z t 1 w b 1[w−ε,w+ε](Ys ) dhY is , (1.5) Lt (Y ) := lim ε→0 2ε 0 where hY i is the predictable quadratic variation process of Y . We define the semimartingale local time for Y at level w ∈ R to be Z 1 t Lw (Y ) = lim 1[w,w+ε](Ys ) dhY is . (1.6) t ε→0 ε 0 It is known (see [2]) that equation (1.4) admits a unique strong solution. Moreover, X is a strong solution to (1.4) if and only if it is a strong solution to     p 1 ρ1 a1 ρ2 a2 1 0 dLt (X) + dLat (X). (1.7) 1− 1− dXt = A(Xt )dBt + 2 ρ2 a2 2 ρ3 a3 The purpose of this paper is to study the parabolic equation (1.2) with initial value u(0, x) = h(x) and solve it as explicitly as possible, where h is a piecewise C N +1 function on R of the form h(x) = h1 (x)1{x≤0} + h2 (x)1{0<x≤a} + h3 (x)1{x>a} , (1.8) with hi being bounded with continuously differentiable up to order N + 1, i = 1, 2, 3, for some integer N ≥ 1. This kind of PDEs appears in the modelization of diffusion phenomena in many areas, for example, in geophysics [9], ecology [4], biology [11] and so on. The non-smoothness of the coefficients reflects the multilayered medium. The solution u(t, x) of (1.2) with initial value h(x) admits a probabilistic representation: u(t, x) = Ex [h(Xt )], t ≥ 0. (1.9) Here the subscript x in Ex means the expectation is taken with respect to the law of the diffusion process X that starts from x at time t = 0. In the first part of this paper, we show that, when √ √ √ √ either ρ1 a1 = ρ2 a2 or ρ2 a2 = ρ3 a3 , we can reduce the diffusion process of (1.3) to skew Brownian motion through transformations by scale functions. Thus in this case, we are able to derive an explicit formula for the transition density function q(t, x, y) of the diffusion process X with respect to the measure m(dx). The kernel q(t, x, y) is the fundamental solution (also called 2 heat kernel) q(t, x, y) of L. Skew Brownian motions are a subclass of diffusion processes X of (1.3) with a1 = a2 = a3 = 1 and ρ2 = ρ3 . It is first introduced by Ito and McKean in 1963, and has since been studied extensively by many authors; see, for example, [1, 2, 13, 15] and the references therein. However we do not know the explicit formula for q(t, x, y) in the general case. In [6], formulas for the fundamental solutions are derived in terms of inverse Fourier transform and Green functions, but they are not explicit. Though we do not have explicit formula for its heat kernel, in the second part of this paper, we are able to derive asymptotic expansion in small time of solutions of (1.2) in the general case with piecewise C N +1 initial data, by making use of the explicit heat kernel obtained under the special case (a2 , ρ2 ) = (a3 , ρ3 ). The small-time asymptotic expansion we get in this paper is explicit and can be directly used in applications. This extends the work of [17] and of [18], where only one discontinuity is studied. In [5] and [10], numerical methods are proposed to study (1.4) with discontinuous coefficients. The rest of the paper is organized as follows. In the section 2, we derive explicit formula for √ √ √ √ the fundamental solution of L in the case when either ρ1 a1 = ρ2 a2 or ρ2 a2 = ρ3 a3 . We then study the small-time asymptotic expansion of the solution of (1.2) in the general case for any positive constants a, ai and ρi . Extensions to the more general piecewise constant coefficient case with more than two discontinuities are mentioned at the end of the paper. 2 Fundamental solution In view of the probabilistic representation (1.9), properties of solution to the heat equation (1.9) can be deduced from the properties of the diffusion process X of (1.7). For notational simplicity, let’s rewrite SDE (1.7) as   ( dYtx = p1{Ytx ≤0} + q1{0<Ytx ≤a} + r1{a<Ytx } dBt + α2 dL0t (Y x ) + β2 dLat (Y x ), (2.1) Y0x = x ∈ R, where p, q, r are positive constants, α, β ∈ (−∞, 1) and B a one-dimensional Brownian motion. The superscript x indicates the process Y x starts from x at t = 0. In the particular case when p = q = r = 1 and β = 0, Ytx is just a Skew Brownian motion of parameter (2 − α)−1 (cf. [15]), and in the case when p = q = r = 1, Ytx is the double-skewed Brownian motion (cf. [12]). The next two theorems are known; see [2, Theorems 2.1 and 2.2] or [8]. We give a proof here not only for reader’s convenience but also some formulas in the proof will be used later to derive the transition density function of Y . Theorem 2.1 For every α < 1 and β < 1, SDE (2.1) has a unique strong solution for every x ∈ R. Proof. Define   x/p s(x) = x/q   (x − a)/r + a/q for x < 0, for x ∈ [0, a], for x > a. Clearly, s is strictly increasing and one-to-one. Let σ denote the inverse of s:   for x < 0, px σ(x) = qx for x ∈ [0, s(a)],   r(x − s(a)) + s(a)q for x > s(a). 3 (2.2) (2.3) Let s′ℓ and σℓ′ denote the left hand derivative of s and σ, respectively; that is,     1/p for x ≤ 0, p for x ≤ 0,  ′ ′ sℓ (x) = 1/q for x ∈ (0, a], and σℓ (x) = q for x ∈ (0, s(a)],     r for x > s(a), 1/r for x > a, The second derivative σ ′′ (dx) = (q − p)δ{0} (dx) + (r − q)δ{s(a)} (dx). Here δ{z} denotes the Dirac measure concentrated at point z. Since both q(α−1) + 1 and r(β−1) + 1 are strictly less than 1, by Theorem 2.1 of [2], the following p q SDE has a unique strong solution Z = {Zt , t ≥ 0} with Z0 = s(x):     1 q(α − 1) 1 r(β − 1) s(a) 0 dZt = dBt + (2.4) + 1 dLt (Z) + + 1 dLt (Z). 2 p 2 q Let Y = σ(Z). We will show that Y is a solution to (2.1). Clearly Y0 = x. Function σ is a piecewise linear function and can be expressed as a difference of two convex functions. By Ito-Tanaka’s formula (see [13, Theorem VI.1.5]), Y is a semimartingale and in fact, for t ≥ 0, Z Z t 1 ′′ Lw σℓ′ (Zs )dZs + Yt = σ(Z0 ) + t (Z)σ (dw) 2 0 Z t  p1{Yt ≤0} + q1{0<Yt ≤a} + r1{Yt >a} dBt = x+ 0 β r s(a) αq 0 Lt (Z) + L (Z). + 2 2 t By Corollary VI.1.9 of Revuz and Yor [13], for t ≥ 0, Z Z 1 t 1 t s(a) Lt (Z) = lim 1[s(a),s(a)+ε) (Zs )dhZis = lim 1[a,a+rε) (Ys )ds ε↓0 ε 0 ε↓0 ε 0 Z 1 t 1 lim = 1 (Ys )dhY is = Lat (Y )/r. r ε↓0 rε 0 [a,a+rε) (2.5) (2.6) A similar calculation shows that L0t (Z) = L0t (Y )/q. Thus we have from (2.5) that Y = σ(Z) is a strong solution for (2.1) with Y0 = x. We now examine pathwise uniqueness. Suppose that Y x is a solution of (2.1). Since s is the difference of two convex functions, we can apply Ito-Tanaka’s formula to get Z Z t 1 x ′′ x x ′ x Lw (2.7) sℓ (Ys )dYs + s(Yt ) = s(x) + t (Y )s (dw). 2 R 0 Let Z = s(Y x ). A similar calculation as above shows that Z satisfies SDE (2.4) with initial value Z0 = s(x). Since by Theorem 2.1 of [2] solutions to (2.4) is unique, hence so is solution Y x to (2.1).  Theorem 2.2 The process Y x is a strong solution of the SDE (2.1), if and only if Y x is a strong solution to   α b 0t (Y x ) + β dL b a (Y x ) dYtx = p1{Ytx ≤0} + q1{0<Ytx ≤a} + r1{a<Ytx } dBt + dL (2.8) 2−α 2−β t with Y0x = x. 4 Proof. This is due to the fact that, for solution Y x of (2.1) (see the proof of Theorem 2.2 in [2]), L0t (Y x ) = 2 b0 x L (Y ), 2−α t Lat (Y x ) = x bw x / {0, a}. and Lw t (Y ) = L (Y ) for w ∈ 2 ba x L (Y ), 2−β t (2.9)  Corollary 2.3 For every p, q > 0, α < 1 and x ∈ R, the following SDE has a unique strong solution for every x ∈ R:   ( dXtx = p1{Xtx ≤0} + q1{0<Xtx } dBt + α2 dL0t (X x ), (2.10) X0x = x ∈ R. Moreover, the transition probability density of the diffusion X x is given by      1{y≤0} 1{y>0} (f (x) − f (y))2 1 X + × exp − p (t, x, y) = √ p q 2t 2πt   p + q(α − 1) (| f (x) | + | f (y) |)2 + sign(y) exp − p − q(α − 1) 2t (2.11) where f (y) = yp 1{y≤0} + yq 1{y>0} . Proof. Taking q = r and β = 0, we have by Theorem 2.1 that for every p, q > 0, α < 1 and x ∈ R,   α (2.12) dXtx = p1{Xtx ≤0} + q1{0<Xtx } dBt + dL0t (X x ), 2 has a unique strong solution with X0x = x. Let s be the function defined by (2.2). We see from (2.4) that Z := s(X x ) satisfies SDE (2.4) with q = r and β = 0; that is,   1 q(α − 1) + 1 dL0t (Z). (2.13) dZt = dBt + 2 p Thus Z is a skew Brownian, whose transition density function is explicitly known. Using Theorem b 0 (Z): 2.2, we can rewrite Z in terms of its symmetric semimartingale local time L dZt = dBt + p + q(α − 1) b 0 dL (Z). p − q(α − 1) t (2.14) Thus by [15] (see also Exercise III.(1.16) on page 82 of [13]), Z has a transition density function pZ (t, x, y) with respect to the Lebesgue measure on R given by      1 (|x| + |y|)2 p + q(α − 1) (x − y)2 Z p (t, x, y) = √ sign(y) exp − + . exp − 2t p − q(α − 1) 2t 2πt The desired formula (2.11) now follows from the fact that pX (t, x, y) = 1 pZ (t, s(x), s(y)) σl′ (s(y))  5 Remark 2.4 (i) When α = 1−(p2 /q 2 ), Corollary 2.3 in particular recovers the formula given in [1, Theorem 5.1] (see the proof of Proposition 5.1 of [1] for a more explicit one) for the transition density function of X x in (2.10). Note that in this case, the more compact formula (2.11) in this paper is equivalent to the form given in [1]. Using the connection to skew Brownian motion established in the proof of Corollary 2.3, one can derive the joint distribution of X x given by (2.8), its semimartingale (respectively, symmetric semimartingale) local time b 0 (X x )) at level 0 and its occupation time on [0, ∞) from that of the L0 (X x ) (respectively, L corresponding skew Brownian motion. The latter can be found in [1, Theorem 1.2]. Same remark applies to the diffusion process Y x satisfying SDE (2.16) below. (ii) The diffusion process X of (2.10) can be called oscillating skew Brownian motion, or skewed oscillating Brownian motion as suggested by Wellner [16]. When α = 1−(p2 /q 2 ), the marginal distribution of X x of (2.10) (characterized by its transition density function pX (t, x, y)) is the Fechner distribution in statistics; see [16]. This corresponds exactly to the case when the density function pX (t, 0, y) is continuous in y. For generator L in (1.1) or its associated SDE (1.7), the condition α = 1 − (p2 /q 2 ) corresponds precisely to the case of ρ2 = ρ3 .  The next theorem allows us to obtain the transition density function for solution Y of (2.1) when β = 1 − rq . Theorem 2.5 Let p, q, r, a > 0 and α < 1. For x ∈ R, let X x denote the unique strong solution of (2.10) with initial value X0x = x. The process Y defined by:   if Y0 = x ≤ a a + qr (Xtx − a)+ − (Xtx − a)− Yt = (2.15) q q  a + r (X r (x−a)+a − a)+ − (X r (x−a)+a − a)− if Y0 = x > a. t t q is the unique strong solution of the stochastic differential equation:   q a x 1 α 1− dLt (Y ) dYtx = p1{Ytx ≤0} + q1{0<Ytx ≤a} + r1{a<Ytx } dBt + dL0t (Y x ) + 2 2 r (2.16) with initial value Y0x = x. Proof. We will present the proof just in the case where x ≤ a. The other case is similar. Consider the bijective function defined by: ϕ(x) = a + (r/q)(x − a)+ − (x − a)− for x ∈ R. Then Y = ϕ(X x ). Note that ϕ(x) = x for x ≤ a, and ϕ(x) > a if and only if x > a. Since ϕ is the difference of two convex functions, applying the Ito-Tanaka formula to Y = ϕ(X x ), we obtain: Z t Z 1 x ′′ ϕ′ℓ (Xsx )dXsx + Yt = Y0 + Lw t (X )ϕ (dx) 2 0 R   Z t  1 r r x − 1 Lat (X x ) 1{Xsx ≤a} + 1{Xsx >a} dXs + = x+ q 2 q 0 Z t  p1{Xsx ≤0} + q1{0<Xsx ≤a} + r1{Xsx >a} dBt = x+ 0 6   α 0 x 1 r + dLt (X ) + − 1 Lat (X x ) 2 2 q Z t  p1{Ysx ≤0} + q1{0<Ysx ≤a} + r1{Ysx >a} dBt = x+ 0   α 0 x 1 r + dLt (X ) + − 1 Lat (X x ). 2 2 q (2.17) By the same reasoning as that for (2.6), we have q L0t (X x ) = L0t (Y ) and Lat (X x ) = Lat (Y ). r (2.18) This together with (2.17) shows that Y is a strong solution to SDE (2.16) with Y0 = x. The uniqueness follows from Theorem 2.1.  Corollary 2.6 When β = 1 − qr , the solution Y to SDE (2.1) has a transition probability density pY (t, x, y) with respect to the Lebesgue measure on R given by ! ( ! 2 1 1 1 1 (s(x) − s(y)) {0<y≤a} {y>a} {y≤0} + + pYt (x, y) = √ × exp − p q r 2t 2πt (2.19) !) (| s(x) | + | s(y) |)2 p + q(α − 1) sign(y) exp − + p − q(α − 1) 2t Proof. This follows immediately from Theorem 2.5 and Corollary 2.3.  √ √ Remark 2.7 SDE (2.1) with β = 1 − qr corresponds exactly to SDE (1.7) with ρ2 a2 = ρ3 a3 . If we use q(t, x, y) to denote the transition density function of Y with respect to its symmetrizing measure m(dx) = ρ(x)dx, then clearly we have q(t, x, y) = pY (t, x, y)/ρ(y). 3 Small-time asymptotic expansion In this section, we will give an explicit asymptotic expansion, in small time, of the function u(x, t) = E(h(Ytx )), in the general case where p, q, r > 0 and α, β ∈ (−∞, 1). Definition 3.1 For every k ∈ N, we denote by erfck the function defined on R by: Z +∞ 2 2 uk e−u du, z ∈ R. erfck (z) = √ π z   Observe that erfck (0) = √1π Γ k+1 , where Γ is the Euler Gamma function. The following 2 elementary properties are straightforward from its definition. Lemma 3.2 For every given integer n ≥ 1, 7 (i) As z → +∞, erfck (z) = o(z −n ).     (ii) As z → −∞: erfck (z) = √1π 1 + (−1)k Γ k+1 + o(|z|−n ). 2 o(|z|−n ) =0 z→+∞ |z|−n Here the notation o(|z|−n ) as z → +∞ (respectively, z → −∞) means that lim o(|z|−n ) = 0). Similar notation will also be used for o(tn ) as t → 0. Note that z→−∞ |z|−n by Corollary 2.3, SDE   ( dStx = q1{Stx ≤a} + r1{a<Stx } dBt + β2 dLat (S x ) (3.1) S0x = x ∈ R (respectively, lim has a unique strong solution for every x ∈ R and it has a transition probability density with respect to the Lebesgue measure on R given by ! ( ! 1{y≤a} 1{y>a} 1 (g(x) − g(y))2 p(t, x, y) = √ × exp − + q r 2t 2πt (3.2) !) q + r(β − 1) (|g(x)| + |g(y)|)2 , + sign(y − a) exp − q − r(β − 1) 2t where g(x) = xq 1{x≤a} + xr 1{x>a} . The following estimation is a key in small-time asymptotic expansion for the solution of heat equation (1.2). It reduces the case to diffusions of the form (2.10) for which we have explicit knowledge about their heat kernels. Lemma 3.3 Define  x  X , solution of sde (2.10) a 2 Ξx = a  S x , solution of sde (3.1) if x > . 2 Then there exist positive constants c1 and c2 such that for every t > 0,     P Ξs 6= Ysx for some s ∈ [0, t] ≤ c1 exp − c2 a2 /t , if x ≤ (3.3) where Y x is the diffusion process of (2.1). a Proof. We will prove the lemma in the case where x ≤ . The other case is similar. For t > 0, let 2 us denote by τt = inf{t ≥ 0; Ξs = a} = inf{t ≥ 0; Xsx = a}, and n o n o At = ω ∈ Ω : Ξs (ω) 6= Ysx (ω) for some s ∈ [0, t] = ω ∈ Ω : Xsx (ω) 6= Ysx (ω) for some s ∈ [0, t] . (3.4) 8 It is clear that, At ⊂ {τt ≤ t} and consequently,     x x P(At ) ≤ P(τt ≤ t) ≤ Px sup |Xs − x| ≥ |a − x| ≤ Px sup |Xs − x| ≥ a/2 . 0≤s≤t 0≤s≤t Here Px is the probability law of X x . Since X is a diffusion whose transition density function q(t, x, y) with respect to the measure m(dx) = ρ(x)dx enjoys the Aronson type Gaussian estimate (1.3), we have by a similar argument as that for [14, Lemma II.1.2] that   2 Px sup |Xsx − x| ≥ a/2 ≤ c1 e−c2 a /t . 0≤s≤t This proves the lemma.  Corollary 3.4 Let Ξx be defined by (3.3). Then for every integer n ≥ 1 and time t > 0 small, i h i h E h(Ytx ) = E h(Ξxt ) + o(tn ). a ; the other case is similar. For all t ≥ 0, 2 i i h h i h E h(Ytx ) = E h(Ytx ); At + E h(Ytx ); Act , Proof. Let us consider the case where x ≤ where At is the set defined in (3.4), and Ac is the complementary of A. Since Xsx (w) = Ysx (w) for all s ∈ [0, t] on Act , we have i i h h i h E h(Ytx ) = E h(Xtx ); Act + E h(Ytx ); At i h i h = E h(Xtx ) + E h(Ytx ) − h(Xtx ); At . The desired result now follows from Cauchy-Schwarz inequality, the fact that h is bounded and Lemma 3.3.  Let h be a piecewise C N +1 function of the form (1.8). By Taylor expansion, for every x ∈ R \ {0, a}, N X h(j) (x) (y − x)j + O(|y − x|N +1 ), h(y) = j! j=1 and for x ∈ {0, a} and i = 1, 2, 3, hi (y) = N (j) X h (x) i j=1 j! (y − x)j + O(|y − x|N +1 ). Here the notation O(|y − x|N +1 ) means that there are constants C, δ > 0 so that the term O(|y − x|N +1 ) is no larger than C|y − x|N +1 for any y ∈ R with |y − x| < δ. Thus we have the following. 9 Proposition 3.5 If x ∈ R \ {0, a}, E [h(Ξxt )] = N  i  h X 1 (j) h (x)E (Ξxt − x)j + O |Ξxt − x|N +1 , j! (3.5) j=0 and if x ∈ {0, a}, E[h(Ξxt )] = N N i X i h h X 1 (j) 1 (j) h1 (x)E (Ξxt − x)j 1{Ξxt ≤0} + h2 (x)E (Ξxt − x)j 1{0<Ξxt ≤a} j! j! j=0 + j=0 N  i  h X 1 (j) h3 (x)E (Ξxt − x)j 1{Ξxt >a} + O |Ξxt − x|N +1 . j! (3.6) j=0 Using the transition probability densities of the diffusions X x and S x given in corollary 2.3 and in (3.2), we can compute the expectations in (3.5)-(3.6). Proposition 3.6 Let k ≥ 1 be an integer. (i) For x < 0, k   i  x  X h  −x  k j/2−1 j/2 k k k k/2−1 k/2 x t erfck √ E (Ξt − x) = (−1) p 2 2 t Aj (x) erfcj √ + , j p 2t p 2t j=0 where   q xk−j (p+q(α−1))pj (−1)k+1 2k−j +2pq j ( −1)k−j Aj (x) = p − q(α − 1) p and   k k! . = j j!(k − j)! (ii) For 0 < x < a, E h (Ξxt − x) k i k/2−1 k k/2 =2 q t k    x   −x  X k j/2−1 j/2 , 2 t Bj (x) erfcj √ + erfck √ j q 2t q 2t j=0 where Bj (x) = p  k−j  xk−j − 2q(α − 1)(−p)j −1 + (p + q(α − 1))q j (−2)k−j . p − q(α − 1) q (iii) For x > a, k   h i a − x X x − a k j/2−1 j/2 k x k/2−1 k k/2 E (Ξt − x) = 2 r t erfck √ 2 t Cj (x) erfcj √ + , j r 2t r 2t j=0 where Cj (x) = q k−j  (x − a)k−j  − 2r(β − 1)(−1)j q j −1 + (q + r(β − 1)) r j 2k−j (−1)k−j . q − r(β − 1) r As of the above Proposition, and by Lemma 3.2, we obtain the following expansion i h a consequence k x of E (Ξt − x) for x ∈ R \ {0, a} and small time t > 0. 10 Corollary 3.7 For any positive integers n > k ≥ 1, x ∈ R \ {0, a}, as t → 0+ , i h E (Ξxt − x)k = 2k/2−1 Dk (x)tk/2 + o(tn ), where 1 Dk (x) = √ (1 + (−1)k )Γ π  k+1 2 h i pk 1{x<0} + q k 1{0<x<a} + r k 1{x>a} . (3.7) Using again the transition probability densities of the diffusions X x and S x , we get Proposition 3.8 For x ∈ {0, a} and integer k ≥ 1,  q(1 − α) (2t)k/2  k + 1     when x = 0 (−1)k pk √ Γ i h p − q(α − 1) 2 π E (Ξxt − x)k 1{Ξxt ≤0} =   r(1 − β) a   when x = a, (−1)k q k (2t)k/2 erfck √  q − r(β − 1) q 2t  i  a  h 1 k + 1 k k k/2 x √ Γ − erfck √ E (Ξt − x) 1{0<Ξxt ≤a} = Ak (x) q (2t) , 2 π q 2t   a  p k k/2 k/2  √  when x = 0 q 2 t erfc k  p − q(α − 1) i h q 2t k x E (Ξt − x) 1{Ξxt >a} = q 2k/2 tk/2  k + 1    when x = a. rk √ Γ  q − r(β − 1) π 2 Here Ak (x) =  p     p − q(α − 1)  r(1 − β)(−1)k    q − r(β − 1) when x = 0 (3.8) when x = a. From Proposition 3.8 and Lemma 3.2 we deduce: Corollary 3.9 For x ∈ {0, a} and integers n > k ≥ 1,   k/2 k/2    q(1 − α) (−1)k pk 2 √t Γ k + 1 i h k x E (Ξt − x) 1{Ξxt ≤0} = p − q(α − 1) 2 π  o(tn ) i h 1 k + 1 E (Ξxt − x)k 10<Ξxt ≤a = Ak (x)q k (2t)k/2 √ Γ + o(tn ), 2 π   when o(tn ) i h k x   k/2 k/2 E (Ξt − x) 1Ξxt >a = q 2 t k+1  when rk √ Γ  q − r(β − 1) 2 π where Ak (x) is defined in (3.8). Combining Proposition 3.5 with Corollaries 3.7 and 3.9, we have Theorem 3.10 For small time and x ∈ R, E [h(Ξxt )] = N X bk (x)tk/2 + O(t(N +1)/2 ), k=0 11 when x = 0 when x = a; x=0 x = a, where  1 ∂ k h(x)  k/2−1   2 D (x) k   k! ∂xk  k    ∂ h1 (0) k+1 ∂ k h2 (0) k 2k/2 k k √ Γ( ) q(1 − α)(−1) p + pq bk (x) = k!(p − q(α − 1)) π 2  ∂xk ∂xk      ∂ k h2 (a) ∂ i h3 (a) k 2k/2 k+1  k k  √ Γ( ) r(1 − β)(−1) q + qr  2 ∂xk ∂xk k!(q − r(β − 1)) π when x ∈ R \ {0, a} when x = 0 when x = a (3.9) with the function Dk (x) given by (3.7). By Theorem 3.10 and Corollary 3.4 we get: Corollary 3.11 For every x ∈ R, when t > 0 is small, E[h(Ytx )] = N X bk (x)tk/2 + O(t(N +1)/2 ), k=0 where the function bk (x) is defined by (3.9). Remark 3.12 (i) When a = 0 and r = q, we recover the expansion obtained in [17]. (ii) Employing the same approach used in this section, one can obtain without any difficulty a similar small-time expansion for the solution of the heat equation (1.2) where A(x) and ρ(x) are piecewise constant with n ≥ 3 points of discontinuity; that is, in case of the following SDE:  n q X αi ai x  dY x = A(Y x )dBt + dLt (Y ), t t 2 i=1   Y0 = x ∈ R, where a1 = 0 < a2 < ... < an , pi ∈ (0, +∞), αi ∈ (−∞, 1) for i = 1, 2, ..., n and A(Ytx ) = p0 1{Ytx ≤a1 } + n−1 X i=1 pi 1{ai <Ytx ≤ai+1 } + pN 1{an <Ytx } . See [2] for the existence and uniqueness of strong solution of this SDE. References [1] T. Appuhamillage, V. Bokil, E. Thomann, E. Waymire and B. Wood, Occupation and local times for skew Brownian motion with applications to dispersion across an interface. Ann. Appl. Probab. 21 (2011), 183-214. Corrections: Ann. Appl. Probab. 21 (2011), 2050-2051. [2] R. F. Bass and Z.-Q. Chen, One-dimensional stochastic differential equations with singular and degenerate coefficients. Sankhya 67 (2005), 19-45. [3] K. Burdzy and Z.-Q. Chen, Local time flow related to skew Brownian motion. Ann. Probab. 29 (2001), 1693-1715. 12 [4] R. Cantrell, C. Cosner, Diffusion models for population dynamics incorporating individual behavior at boundaries: Applications to refuge design. Theor. Population Biology, 55, (1999), 198-207. [5] P. Étoré, On random walk simulation of one-dimonsional diffusion processes with discontinuous coefficients. Electron. J. Probab. 11 (2006), 249-275. [6] B. Gaveau, M. Okada, T. Okada, Second order differential operators and Dirichlet integrals with singular coefficients. Tohoku Math. J. 39 (1987), 465-504. [7] P. Gyrya and L. Saloff-Coste. Neumann and Dirichlet Heat Kernels in Inner Uniform Domains. Paris: Socit mathmatique de France, 2011. [8] J. G. Le Gall, One-dimensional stochastic differential equations involving the local times of the unknown process. In Stochastic Analysis and Applications (Swansea, 1983). Lecture Notes in Math. 1095 (1984), 5182. Springer, Berlin. [9] A. Lejay, Monte Carlo methods for fissured porous media: a gridless approach. Monte Carlo Methods Appl. 10, (2004), 385-392. [10] A. Lejay and M. Martinez, A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients. Annals Appl. Probab. 16 (2006), 107-139. [11] S. Nicas, Some results on spectral theory over networks, applied to nerve impulse transmission. Orthoginal Polynomials and Applications (Bar-le-Duc, 1984). Lect. notes Math. 1171. Springer, 532-541. [12] J. M. Ramirez, Multi-skewed Brownian motion and diffusion in layered Media. Proc. Amri. Math. Soc. 139 (2011), 3739-3752. [13] D. Revuz and M. Yor, Continuous Martingales and Brownian Motion. Springer-Verlag, Berlin, 1991. [14] D. W. Stroock, Diffusion semigroups corresponding to uniformly elliptic divergence form operators. Séminaire de Probabilités, XXII, 316-347, Lecture Notes in Math., 1321, Springer, Berlin, 1988. [15] J. Walsh, Diffusion with discontinuous local time. Astérisques 52-53 (1978), 37-45. [16] J. Wellner, Fechner’s distribution and connections to skew Brownian motion. Preprint, 2013. [17] M. Zili, Développement asymptotique en temps petits de la solution d’une équation aux dérivées partielles de type parabolique généralisée au sens des distributions-mesures. Note des Comptes Rendues de l’Académie des Sciences de Paris, t. 321, Série I, p. 1049-1052, 1995. [18] M. Zili, Construction d’une solution fondamentale d’une équation aux dérivées partielles à coefficients constants par morceaux. Bull. Sci. Math. 123 (1999), 115-155. Zhen-Qing Chen Department of Mathematics, University of Washington, Seattle, WA 98195, USA E-mail: zqchen@uw.edu Mounir Zili Department of Mathematics, Faculty of sciences of Monastir, Tunisia E-mail: Mounir.Zili@fsm.rnu.tn 13