Qualitative Properties of a Numerical Scheme for the Heat Equation Liviu I. Ignat Departamento de Matem´ aticas, Facultad de Ciencias, Universidad Aut´ onoma de Madrid, 28049 Madrid, Spain [email protected] Summary. In this paper we consider a classical finite difference approximation of the heat equation. We study the long time behaviour of the solutions of the considered scheme and various questions related to the fundamental solutions. Finally we obtain the first term in the asymptotic expansion of the solutions.

1 Introduction The main goal of this paper is the study of the long time behaviour of classical finite difference approximations of the heat equation. Let us consider the linear heat equation on the whole space ( ut − ∆u = 0 in Rd × (0, ∞), u(0, x) = ϕ(x) in Rd . By means of Fourier’s transform, solutions can be represented as the convolutions between the fundamental solutions and the initial data: u(t) = G(t, ·) ∗ ϕ, where

|x|2 1 1 − 4t e = G(t, x) = (2π)d (4πt)−d/2

Z 2

eix·ξ e−|ξ| t dξ. Rd

The smoothing effect of the fundamental solutions G(t, x) yields to the following behaviour of the solution (cf. [3], Ch. 3, p. 44): ku(t)kLp (Rd ) ≤ C(p, q) t−d/2 (1/q−1/p) kϕkLq (Rd ) ,

t > 0, p ≥ q.

(1)

A finer analysis is given in [4], where the authors consider initial data which decay polynomially at infinity. Duoandikoetxea & Zuazua [4] study how the mass of the solution is distributed as t → ∞. They prove the existence of a

2

Liviu I. Ignat

positive constant c = c(p, q, d) such that for any q and p satisfying 1 ≤ q < d/(d − 1), d ≥ 2 (1 ≤ q < ∞ for d = 1), q ≤ p < ∞, ° ° µZ ¶ ° ° d 1 1 1 °u(t, ·) − ° ϕ(x)dx G(t, ·) ≤ c t− 2 − 2 ( q − p ) k|x|ϕkLq (Rd ) (2) ° ° Rd

Lp (Rd )

holds for all t > 0 and ϕ ∈ L1 (Rd ) with |x|ϕ(x) ∈ Lq (Rd ). Let us consider the classical finite-difference scheme:   duh   = ∆h uh , t > 0, dt    uh (0) = ϕh .

(3)

Here uh stands for the infinite unknown vector {uhj }j∈Zd , uhj (t) being the approximation of the solution u at the node xj = jh, and ∆h is the classical second order finite difference approximation of ∆: d 1 X h (uj+ek + uhj−ek − 2uhj ). h2

(∆h uh )j =

k=1

This scheme is widely used and satisfies the classical properties of consistency and stability which imply L2 -convergence (cf. [7], Ch. 13, p. 292). It is interesting to know whenever the properties of the continuous problem are preserved by the numerical scheme. In the following we are concerned with the spatial shape of the discrete solution for large times. To do that we introduce the spaces lp (hZd ): o n X |uj |p < ∞ lp (hZd ) = {uj }j∈Zd : kukplp (hZd ) = hd j∈Zd

and study the behaviour of lp (hZd )-norms of the solutions as t → ∞. The main tool in our analysis is the semi-discrete Fourier transform (SDFT): h π π id X e−ij·ξh uj , ξ ∈ − , u b(ξ) = hd h h d j∈Z

and its inverse uj =

1 (2π)d

Z u b(ξ)eij·ξh dξ, j ∈ Zd . [−π/h,π/h]d

We refer to [5] and [10] for a survey on this subject. By means of SDFT we compute the solutions of equation (3) in a similar way as in the continuous case, writing them as a convolution of a fundamental solution Ktd,h and the initial datum. This allows us to obtain decay rates of the solution in different lq − lp norms analogous to (1). All the estimates are uniform with respect to the step size, h. This proves a kind of lq − lp stability of our scheme:

Qualitative Properties of a Numerical Scheme for the Heat Equation

3

Theorem 1. Let 1 ≤ q ≤ p ≤ ∞. Then there exists a positive constant c(p, q, d) such that kuh (t)klp (hZd ) ≤ c(p, q, d) t−d/2 (1/q−1/p) kϕh klq (hZd ) for all t > 0, uniformly in h > 0. A similar approach in the case of the transport equation has been studied by Brenner and Thom´ee [2] and Trefethen [11]. They introduce a finite difference approximation and give conditions which guarantee the lp -stability of the scheme. Next we prove that the fundamental solutions Ktd,h of equation (3) are related to the modified Bessel function Iν (x): µ (Ktd,h )j =

exp(− h2t2 ) πh

¶d Y d

µ Ijk

k=1

2t h2

¶ , j = (j1 , j2 , . . . , jd ) ∈ Zd .

(4)

This property proves the positivity and various properties regarding the monotonicity of the discrete kernel Ktd,h . Finally, we consider the weighted space l1 (hZd , |x|) and obtain the first term in the asymptotic expansion of the discrete solution. The weighted spaces lp (hZd , |x|), 1 ≤ p < ∞ are defined as follows: o n X |uj |p |jh|p < ∞ . lp (hZd , |x|) = {uj }j∈Zd : kukplp (hZd ,|x|) = hd j∈Zd

The following theorem gives us the first term of the asymptotic expansion of the solution uh : Theorem 2. Let p ≥ 1. Then there exists a positive constant c(p, d) such that ° °   ° ° X ° h ° d,h ° h °u (t) − h ϕ j Kt ° ≤ c(p, d) t−1/2−d/2 (1−1/p) kϕh kl1 (hZd ,|x|) ° ° ° d j∈Z p d l (hZ )

for all ϕh ∈ l1 (hZd , |x|) and t > 0, uniformly in h > 0. This shows that for t large enough the solution behaves as the fundamental solution. In contrast with (2) our result is valid only for the initial data in the weighted space l1 (hZ, |x|). The extension of this result to general initial data, i.e. in lq (hZ, |x|), 1 < q < p, remains an open problem. In [6] we consider the first k ≥ 1 terms of the asymptotic expansion of the discrete solution and obtain a similar result.

2 Proof of the Results By means of SDFT we obtain that u bh satisfies the following ODE:

4

Liviu I. Ignat

µ ¶ d h π π id db uh 4 X 2 ξk h (t, ξ) = − 2 sin u bh (t, ξ), t > 0, ξ ∈ − , . dt h 2 h h k=1

In the Fourier space, the solution u bh reads h π π id u bh (t, ξ) = e−tph (ξ) ϕ bh (ξ), ξ ∈ − , , h h where the function ph : [−π/h, π/h]d → R is given by µ ¶ d 4 X 2 ξk h ph (ξ) = 2 sin . h 2

(5)

k=1

The solution of equation (3) is given by a discrete convolution between the fundamental solution Ktd,h and the initial datum: uh (t) = Ktd,h ∗ ϕh . The inverse SDFT of the function e−tph (ξ) gives us the following representation of the fundamental solution Ktd,h : Z 1 (Ktd,h )j = e−tph (ξ) eij·ξh dξ, j ∈ Zd . (2π)d [−π/h,π/h]d We point out that for any j = (j1 , j2 , . . . , jd ) ∈ Zd the kernel Ktd,h can be written as the product of one-dimensional kernels Kt1,h : (Ktd,h )j

=

d Y

(Kt1,h )jk .

(6)

k=1

A simple change of variables in the explicit formula of Kt1,h relates it with the modified Bessel functions: µ ¶ exp(− h2t2 ) 2t (Kth )j = Ij , j ∈ Z. πh h2 Separation of variables formula (6) proves (4). We recall that the modified Bessel’s function Iν (x) is positive for any positive x. Also for a fixed x, the map ν → Iν (x) is even and decreasing on [0, ∞) (cf. [8], Ch. II, p. 60). These properties prove that the kernel Ktd,h has the following properties: Theorem 3. Let t > 0 and h > 0. Then i) For any j = (j1 , j2 , . . . , jd ) ∈ Zd µ (Ktd,h )j =

exp(− h2t2 ) πh

¶d Y d k=1

µ Ijk

2t h2

¶ .

Qualitative Properties of a Numerical Scheme for the Heat Equation

5

ii) For any j ∈ Zd , the kernel (Ktd,h )j is positive. iii) The map j ∈ Z 7→ (Kt1,h )j is increasing for j ≤ 0 and decreasing for j ≥ 0. iv) For any a = (a1 , a2 , . . . , ad ) ∈ Zd and b = (b1 , b2 , . . . , bd ) ∈ Zd satisfying |a1 | ≤ |b1 |, |a2 | ≤ |b2 |, . . . , |ad | ≤ |bd |, the following holds

(Ktd,h )b ≤ (Ktd,h )a .

The long time behaviour of the kernel Ktd,h is similar to the one of its continuous counterpart. Theorem 4. Let p ∈ [1, ∞]. Then there exists a positive constant c(p, d) such that kKtd,h klp (hZd ) ≤ c(p, d) t−d/2 (1−1/p) (7) holds for all positive times t, uniformly on h > 0. Once Theorem 4 is proved, Young’s inequality provides the decay rates of the solutions of equation (3) as stated in Theorem 1. d,1 Proof (of Theorem 4). A scaling argument shows that (Ktd,h )j = (Kt/h 2 )j , reducing the proof to the case h = 1. In the sequel we consider the band limited interpolator of the sequence Ktd,1 (cf. [12], Ch. I, p. 13): Z 1 K∗d (t, x) = eix·ξ e−tp1 (ξ) dξ. (8) (2π)d [−π,π]d

In [9] the authors prove the existence of a positive constant A such that for any function f with its Fourier transform supported in the cube [−π, π]d the following holds: Z X p d |f (j)| ≤ A |f (x)|p dx, p ≥ 1. (9) j∈Zd

Rd

This reduces (7) to similar estimates on the Lp (Rd )-norm of K∗d . The interpolator K∗d satisfies kDα K∗d (t, ·)kLp (R) ≤ c(α, p, d) t−|α|/2−d/2 (1−1/p)

(10)

for any multiindex α = (α1 , . . . , αd ) and 1 ≤ p ≤ ∞. Using (5) and (8), we reduce (10) to the one dimensional case. We consider the cases p = 1 and p = ∞. The general case, 1 < p < ∞, follows by the H¨older inequality. The case p = ∞ easily follows by the rough estimate: µ ¶ Z π 1 2 ξ α α 1 |ξ| exp −4t sin dξ ≤ c(α) t−(α+1)/2 . kD K∗ (t, ·)kL∞ (Rd ) ≤ 2π −π 2

6

Liviu I. Ignat

Finally, we apply Carlson-Beurling’s inequality (cf. [1] and [2]): kb akL1 (R) ≤ (2kakL2 (R) ka0 kL2 (R) )1/2 to the function a(ξ) = |ξ|α exp(−4t sin2 ξ/2). We obtain the existence of a positive constant C such that for all t > 0, kK∗1 (t)kL1 (R) ≤ C. This proves Theorem 4.

u t

Now we sketch the proof of Theorem 2. Proof (of Theorem 2). First, a scaling argument reduces the proof to the case h = 1. We consider the cases p = 1 and p = ∞, the other cases follow by interpolation. The solution u1 (t) of equation (3) is given by: X u1j (t) = (Ktd,1 ∗ ϕ1 )j = (Ktd,1 )j−n ϕ1n . n∈Zd

Let us introduce the sequence {aj (t)}j∈Zd as follows ³ ´ X X aj (t) = u1 (t) − Ktd,1 ϕ1n = u1j (t) − (Ktd,1 )j ϕ1n =

(Ktd,1 )j−n ϕ1n

n∈Zd

=

j

n∈Zd

X X

ϕ1n

³

−

(Ktd,1 )j

X

n∈Zd

ϕ1n

n∈Zd

(Ktd,1 )j−n

−

(Ktd,1 )j

´

.

n∈Zd

In the sequel we denote by c a constant that may change from one line to another. It remains to prove that sup |aj (t)| ≤ c t−(d+1)/2 kϕ1 kl1 (Zd ,|x|)

j∈Zd

and

X

|aj (t)| ≤ c t−1/2 kϕ1 kl1 (Zd ,|x|) .

j∈Zd

The Taylor formula applied to the function K∗d gives us Z K∗d (t, j − n) − K∗d (t, j) =

1

X

0 |α|=1

Dα K∗d (t, j − sn)(−n)α ds.

As a consequence, for any j ∈ Zd the sequence aj (t) satisfies

(11)

Qualitative Properties of a Numerical Scheme for the Heat Equation

|aj (t)| ≤

X

|ϕ1n |

n∈Zd

≤c

X Z |α|=1

X

|ϕ1n ||n|

n∈Zd

=c

X

1

|nα ||Dα K∗d (t, j − sn)|ds

0

X Z

|α|=1

X

|ϕ1n ||n|

n∈Zd

7

1 0

|Dα K∗d (t, j − sn)|ds

bα j,n (t).

(12)

|α|=1

To prove inequality (11), which corresponds to p = ∞, it is sufficient to show that −(d+1)/2 bα j,n (t) ≤ c t for all indices α with |α| = 1. Inequality (10) shows that α d −|α|/2−d/2 bα = c t−(d+1)/2 . j,n (t) ≤ kD K∗ (t)kL∞ (Rd ) ≤ c t

Now let us consider the case p = 1. We sum on j ∈ Zd in inequality (12) and obtain: X X X X |aj (t)| ≤ |ϕ1n ||n| bα j,n (t) j∈Zd

j∈Zd n∈Zd

X

=

|ϕ1n ||n|

|α|=1

X X

bα j,n (t).

|α|=1 j∈Zd

n∈Zd

It remains to prove that X

−1/2 bα j,n (t) ≤ c t

(13)

j∈Zd

for all n ∈ Zd and for any multiindex α with |α| = 1. Using the separation of variables, we get for all j = (j1 , . . . , jd ) ∈ Zd and n = (n1 , . . . , nd ) ∈ Zd , Z bα j,n (t)

=

1

d Y

0 k=1

|Dα K∗1 (t, jk − snk )|ds

and hence, X j∈Zd

Z bα j,n (t) =

1

d Y

0 k=1

≤ sup s∈R

d Y k=1



 X



jk ∈Z

 

X

|Dαk K∗1 (t, jk − snk )| ds  |Dαk K∗1 (t, jk − s)| .

jk ∈Z

We prove that each term in the last product is dominated by t−αk /2 and consequently the product will be bounded by t−|α|/2 . Applying (9) to the function K∗1 (t, · − s), each of the above sum satisfies

8

Liviu I. Ignat

X

Z |Dαk K∗1 (t, jk − s)| ≤ c

jk ∈Z

R

|Dαk K∗1 (t, x − s)|dx

Z =c

R

|Dαk K∗1,1 (t, x)|dx ≤ c t−|αk |/2 .

This proves inequality (13) and finishes the proof of Theorem 2. t u

Acknowledgements. The author wishes to thank the guidance of his Ph.D advisor Enrique Zuazua. This work has been supported by the doctoral fellowship AP2003-2299 of MEC (Spain) and the grants MTM2005-00714 of the MEC (Spain), 80/2005 of CNCSIS (Romania), “Smart Systems” of the European Union.

References 1. Beurling A.: Sur les int´egrales de Fourier absolument convergentes et leur application ` a une transformation fonctionnelle. In 9. Congr. des Math. Scand., 345–366 (1939). 2. Brenner P., Thom´ee V.: Stability and convergence rates in Lp for certain difference schemes. Math. Scand., 27, 5–23, (1970). 3. Cazenave T., Haraux A.: An introduction to semilinear evolution equations. Oxford Lecture Series in Mathematics and its Applications. 13. Oxford: Clarendon Press. xiv, (1998). 4. Duoandikoetxea J., Zuazua E.: Moments, masses de Dirac et d´ecomposition de fonctions. (Moments, Dirac deltas and expansion of functions). C. R. Acad. Sci. Paris S´er. I Math., 315, (6), 693–698, (1992). 5. Henrici P.: Applied and computational complex analysis. Volume III: Discrete Fourier analysis, Cauchy integrals, construction of conformal maps, univalent functions. Wiley Classics Library. New York, Wiley, (1993). 6. Ignat L.I.: Ph.D. thesis. Universidad Aut´ onoma de Madrid. In preparation. 7. Iserles A.: A first course in the numerical analysis of differential equations. Cambridge Texts in Applied Mathematics. Cambridge Univ. Press, (1995). 8. Olver F.W.J.: Introduction to asymptotics and special functions. Academic Press, (1974). 9. Plancherel M., P´ olya G.: Fonctions enti`eres et int´egrales de Fourier multiples. II. Comment. Math. Helv., 10, 110–163, (1937). 10. Trefethen L.N.: Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations. http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen /pdetext.html. 11. Trefethen L.N.: On lp -instability and oscillation at discontinuities in finite difference schemes. Advances in Computer Methods for Partial Differential Equations, V, 329–331, (1984). 12. Vichnevetsky R., Bowles J.B.: Fourier analysis of numerical approximations of hyperbolic equations. SIAM Studies in Applied Mathematics, 5. SIAM, (1982).