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).

## Qualitative Properties of a Numerical Scheme for the ...

Let us consider the linear heat equation on the whole space. { ut â âu =0 in Rd Ã (0, ..... First, a scaling argument reduces the proof to the case h = 1. We consider ...

#### Recommend Documents

Dispersive properties of a viscous numerical scheme ...
Dans cette Note on analyse les propriÃ©tÃ©s dispersives de quelques schÃ©mas semi-discrets en .... First, using that ph(Î¾) changes convexity at the point Ï/2h, we choose as initial .... of the initial data Ï â L2(R) such that EÏh â Ï weakly

Qualitative properties of generalized principal ...
[3] H. Berestycki, L. Nirenberg, S.R.S. Varadhan, The principal eigenvalue and maxi- mum principle for second-order elliptic operators in general domains, ...

1 Dispersive Properties of Numerical Schemes for ... - CiteSeerX
Keel, M. and Tao, T., (1998). Endpoint Strichartz estimates, Am. J. Math., ... Methods for Or- dinary and Partial Differential Equations, http://web.comlab.ox.ac.uk.

1 Dispersive Properties of Numerical Schemes for ... - CiteSeerX
alternate schemes preserve the dispersion properties of the continuous model. ... the fact that classical energy methods fail, using these dispersion prop- erties, the numerical solutions of the semi-discrete nonlinear problems are proved to ...

A numerical method for the computation of the ...
Considering the equation (1) in its integral form and using the condition (11) we obtain that sup ..... Stud, 17 Princeton University Press, Princeton, NJ, 1949. ... [6] T. Barker, R. Bowles and W. Williams, Development and application of a.

I. The numerical scheme and generalized solitary waves - UoA Scholar
Jun 13, 2006 - aDepartment of Mathematics, Statistics and Computer Science, University of Illinois at ... After suitable testing of the numerical scheme, it is used to examine the .... h of periodic smooth splines of order r (degree r â 1) on.

Numerical optimization of a grating coupler for the ...
tivity of Ag was included to take into account resistive losses in Ag. The input beam was chosen to be of a Gauss- ian profile with a FWHM diameter of 1 m and ...

A comparison of numerical methods for solving the ...
Jun 12, 2007 - solution. We may conclude that the FDS scheme is second-order accurate, but .... âa2 â cos bt + 2 arctan[Î³(x, t)] + bsin bt + 2 arctan[Î³(x, t)]. â ln.

A Methodology for the Construction of Scheme - IJEECS
Internet QoS has actually shown amplified average ... KGB's Internet-2 overlay network to measure the ex- ... sembled using AT&T System V's compiler built on.

A Methodology for the Construction of Scheme - IJEECS
cation of fiber-optic cables and semaphores. Simi- .... sembled using AT&T System V's compiler built on .... [23] K. J. Anderson, âAn analysis of e-business with.

Compound matrices: properties, numerical issues and ...
May 19, 2008. Abstract. This paper studies the possibility to calculate efficiently compounds of real matrices which have a special form or structure. The usefulness of such an effort lies in the fact that the computation of compound matrices, which

A Numerical Study of the Sensitivity of Cloudy-Scene ... - Sites
Jun 28, 1996 - Our civilization has been the single most important force in recent ... The solar energy absorbed and reflected by the earth occurs .... The design strategy of the experiment incorporated both old and new technology concepts.

numerical evaluation of the vibroacoustic behavior of a ...
Oct 24, 2011 - The analytical solution of rectangular cavity (Blevins, 1979; Kinsler et al., .... The speaker box was placed in contact with the side of the cavity ...

A New Algorithm for Finding Numerical Solutions of ... - IEEE Xplore
optimal control problem is the viscosity solution of its associated Hamilton-Jacobi-Bellman equation. An example that the closed form solutions of optimal ...

A new algorithm for finding numerical solutions of ...
Mar 2, 2009 - Department of Mathematics, Beijing Institute of Technology, Beijing .... Given a running cost L(t, y, u) and a terminal cost Ï(y), the optimal ...

numerical approximation of the boundary control for the ...
Oct 11, 2006 - a uniform time T > 0 for the control of all solutions is required. This time T depends on the size of the domain and the velocity of propagation of waves. In general, any semi-discrete dynamics generates spurious high-frequency oscilla

A Quality of Service Routing Scheme for Packet ...
Abstract. Quality of Service (QoS) guarantees must be supported in a network that intends to carry real-time multimedia traffic effectively. A key problem in providing. QoS guarantees is routing which consists of finding a path in a network that sati

Performance evaluation of a reservation random access scheme for ...
We compute the steady state distribution of the Markov chain. This result is used to ... This work is supported by a University of California MICRO and Pacific-Bell ...

Security of A Multisignature Scheme for Specified ...
Oct 9, 2007 - r Â· gw mod p. If it holds, then the verifier accepts the signature is valid; Rejects, otherwise. 3 Security of Zhang et al.'s Multisignature Scheme.

A Reference Discretization Strategy for the Numerical
the relations held between physical quantities within each theory. Let us call ... It must be said that there are still issues that wait to be clarified in relation to.

Modeling of optical properties for a polydispersion/gas ...
Predictive methods for solving the radiative heat transfer equation require the ... All the other commonly used models make recourse to the flux method, that ...

numerical study of the behaviour for elastic- viscoplastic ...
Abstract : The variation of stress during creep convergence of a horizontal circular galleries excavated in rock salt is studied. Examples are given for rock salt by N. Cristescu ([1], [2]). A non-associated elasto-viscoplastic constitutive equation

NUMERICAL DISPERSIVE SCHEMES FOR THE ...
To recover the dispersive properties of the solutions at the discrete level, we ... nonlinear problems with L2-initial data, without additional regularity hypotheses. ... Project CIT-370200-2005-10 in the PROFIT program and the SIMUMAT project ...