domain: phys math nlin info

text:
Denys Dutykh
email: denys.dutykh@univ-savoie.fr

Carl David

Tolmé Runge

Claude Louis

Marie Henri Navier

‡ George

Gabriel Stokes

A brief introduction to pseudo-spectral methods: application to diffusion problems

Keywords: pseudo-spectral methods, approximation theory, diffusion, parabolic equations MSC: [2010] 65M70, 65N35 (primary), 80M22, 76M22 (secondary) PACS: [2010] 47.11.Kb (primary), 44.35.+c (secondary) Key words and phrases. pseudo-spectral methods, approximation theory, diffusion, parabolic equations

The topic of these notes could be easily expanded into a full one-semester course. Nevertheless we shall try to give some flavour along with theoretical bases of spectral and pseudo-spectral methods. The main focus is made on Fourier-type discretizations, even if some indications on how to handle non-periodic problems via Tchebyshev and Legendre approaches are made as well. The applications presented here are diffusiontype problems in accordance with the topics of the PhD school.

General introduction

A numerical simulation is like sex.

If it is good, then it is great. If it is bad, then it is still better than nothing.

-Dr. D (paraphrasing Dr. Z) People think they do not understand Mathematics, but it is all about how you explain it to them.

-I. M. Gelfand

This document represents a collection of brief notes of lectures delivered by the Author at the International PhD school "Numerical Methods for Diffusion Phenomena in Building Physics: Theory and Practice" which took place in Pontifical Catholic University of Parana (PUCPR, Curitiba, Brazil) in April, 2016. The Author of the present document is grateful to the Organizing Committee of this school (in particular to Professors Nathan Mendes and Marx Chhay) for giving him an opportunity to lecture there. The present document should be considered as a supplementary material to the Lectures delivered by the Author at this PhD school. The exposition below is partially based on [START_REF] Fornberg | A practical guide to pseudospectral methods[END_REF] and some other references mentioned in the manuscript.

This lecture is organized as follows. First, we present some theoretical bases behind spectral discretizations in Section 2. An application to a problem stemming from the building physics is given in Section 4. Finally, we give some indications for the further reading in Section 5. This document contains also a certain number of Appendices directly or indirectly related to spectral methods. For instance, in Appendix A we give some useful identities about Tchebyshev polynomials and in Appendix B we give some flavour of Trefftz methods, which are essentially forgotten nowadays. We decided to include also some history of the diffusion in Sciences in general and in Physics in particular. You will find this information in Appendix C. Finally, we prepared also an Appendix D devoted to the Monte-Carlo methods to simulate numerically diffusion processes. The Author admits that it is not related to spectral methods, but nevertheless we decided to include it since these methods remain essentially unknown in the community of numerical methods for Partial Differential Equations (PDEs).

Introduction to spectral methods

Consider for simplicity a 1D compact * domain, e.g. U = [-1, 1] an Ordinary or Partial Differential Equation (ODE or PDE) on it

L u = g , x ∈ U (2.1)
where u(x, t) (or just u(x) in the ODE case) is the solution † which satisfies some additional (boundary) conditions at x = ±1 depending on the (linear or nonlinear) operator L = L (∂ t , ∂ x , ∂ xx , . . .). Function g(x, t) is a source term which does not depend on the solution u. For example, if (2.1) is the classical heat equation in the free space (i.e. without heat sources), then L ≡ ∂ t -ν ∂ xx , g ≡ 0 , where ν ∈ R is the diffusion coefficient. If the solution u(x) is steady (i.e. time-independent), we deal with an ODE. In any case, in the present document we focus only on the discretization in space. The time discretization is discussed deeper in the lecture of Dr. Marx Chhay.

A brief reminder of some most useful numerical techniques for ODEs is given in Appendix F.

The idea behind a spectral method is to approximate a solution u(x, t) by a finite sum

u(x, t) ≈ u n (x, t) = n k = 0 v k (t) φ k (x), (2.2)
where {φ k (x)} ∞ k=0 is the set of basis functions. In ODE case all v k (t) ≡ const. The main question which arises is how to choose the basis functions? Once the choice of {φ k (x)} ∞ k = 0 is made, the second question appears: how to determine the expansion coefficients v k (t)?

Here we shall implicitly assume that the function u(x, t) is smooth. Only in this case the full potential of spectral methods can be exploited ‡ . The concept of smoothness in Mathematics is ambiguous since there are various classes of smooth functions:

C p (U) ⊇ C ∞ (U) ⊇ A ∞ (D), (p ∈ Z + , U ⊆ D ⊆ C).
The last sequence of inclusions needs perhaps some explanations:

C p (U): the class of functions f : U → R having at least p 1 continuous derivatives. * A domain U ⊆ R d , d 1 is compact if it is bounded and closed. For a more general definition of compactness we refer to any course in General Topology.

† For simplicity we consider in this Section scalar equations only. The generalizations for systems is straightforward, since one can apply the same discretization for every individual component of the solution.

‡ We have to say that pseudo-spectral methods can be applied to problems featuring eventually discontinuous solutions as well (e.g. hyperbolic conservation laws). However, it is out of scope of the present lecture devoted rather to parabolic problems. As a side remark, the Author would add that the advantage of pseudo-spectral methods for nonlinear hyperbolic equations is not so clear comparing to modern high-resolution shock-capturing schemes [START_REF] Van Leer | Upwind and High-Resolution Methods for Compressible Flow: From Donor Cell to Residual-Distribution Schemes[END_REF].

C ∞ (U): the class of functions f : U → R having infinitely many continuous derivatives.

A ∞ (D): the class of functions f : D → C analytical (holomorphic) in a domain D ⊆ C containing the segment U in its interior.

For example, the Runge * function f (x) = 1 1 + 25 x 2 is in C ∞ [-1, 1] , but not analytical in the complex plain A ∞ (C).

Choice of the basis

A successful expansion basis meets the following requirements:

(1) [Convergence] The approximations u n (x, t) should converge rapidly to u(x, t) as n → ∞ (2) [Differentiation] Given coefficients {v k (t)} n k=0 , it should be easy to determine another set of coefficients † {v ′ k (t)} n k=0 such that

∂u n ∂x = n k = 0 v k (t) dφ k (x) dx n k = 0 v ′ k (t) φ k (x) .
(3) [Transformation] The computation of expansion coefficients {v k } n k=0 from function values u(x i , t) n i=0 and the reconstruction of solution values in nodes from the set of coefficients {v k } n k=0 should be easy, i.e. the conversion between two data sets is algorithmically efficient u(x i , t)

n i=0
⇆ {v k } n k=0 .

Periodic problems

For periodic problems it is straightforward to propose a basis which satisfies the requirements (1)-(3) above. It consists of trigonometric polynomials:

u n (x, t) = a 0 (t) + n k = 1 a k (t) cos(kπx) + b k (t) sin(kπx) .
(2.

3)

The first two points are explained in elementary courses of analysis, while the requirement (3) is possible thanks to the invention of the Fast Fourier Transform (FFT) algorithm first by Gauß and later by Cooley & Tukey in 1965 [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF].

Remark 1. The use of trigonometric bases such as (2.3) or {e i kπx } k∈Z for heat conduction problems has been initiated by J. Fourier ‡ (1822) [START_REF] Fourier | Théorie analytique de la chaleur[END_REF] even if Fourier series have been known well before Fourier. There is the so-called Arnold's * principle which states that in Mathematics nothing is named after its true inventor. The question the Author would like to rise is whether Arnold's principle is applicable to it-self ?

Remark 2. The term 'Spectral methods' can be now explained. It comes from the fact that the solution u(x, t) is expanded into a series of orthogonal eigenfunctions of some linear operator L (with partial or ordinary derivatives). In this way, the numerical solution is related to its spectrum, thus justifying the name 'spectral methods'. For example, if we take

the Laplace operator L = -∇ 2 ≡ -d j = 1 ∂ 2 ∂x 2 j on the periodic domain [ 0, 2 π) d , its
spectrum consists of Fourier modes:

-∇ 2 e -i k•x = | k | 2 e -i k•x .
In this way we obtain naturally the Fourier analysis and Fourier-type pseudo-spectral methods.

Non-periodic problems

The trigonometric basis (2.3) fails to work for general non-periodic problems essentially because of the failure of requirement [START_REF] Benjamin | Combinatorially composing Chebyshev polynomials[END_REF]. Indeed, the artificial discontinuities arising after the periodisation (see Figure 1) make the Fourier coefficients v n decay as O(n -1 ) when n → ∞.

The analytic series (or Taylor-like expansions) represent another interesting alternative:

u n (x, t) = n k = 1 v k (t) x k .
(2.4)

The problem is that Taylor-type expansions (2.4) work well only for a very limited class of functions. Namely, they satisfy the requirement (1) only for functions analytical † in the unit disc D 1 (0) ⊆ C. For instance, the celebrated Runge's function f (x) = 1 1 + 25x 2 fails to satisfy this condition because of two imaginary poles located at z = ± i 5 . The successful examples of polynomial bases are given by orthogonal polynomials, which arise naturally in various contexts: Numerical integration: Optimal numerical integration formulas achieve a high accuracy by using zeros of orthogonal polynomials as nodes.

Sturm-Liouville problem: Jacobi * polynomials arise as eigenfunctions to singular Sturm † -Liouville ‡ problems.

Approximation in L 2 : Truncated expansions in Legendre § polynomials are optimal approximants in the L 2 -norm.

Approximation in L ∞ : Truncated expansions in Tchebyshev ¶ polynomials are optimal approximants in the L ∞ -norm.

Each of the topics above deserves a separate course to be covered. Here we content just to provide this information as facts, which can be deepened later, if necessary. We would just like to quote Bernardi & Maday (1991):

We do think that the corner stone of collocation techniques is the choice of the collocation nodes [ . . . ] in spectral methods these are always built from the nodes of a Gauß quadrature formula. Consequently, extrema (and zeros) of Tchebyshev (and some other orthogonal) polynomials play a very important rôle in the Numerical Analysis (NA). Tchebyshev nodes are given explicitly by

x k = -cos πk N ∈ [-1, 1] , k = 0, 1, 2, . . . , N . (2.5)
There is a result saying that using nodes (2.5) as interpolation points gives an interpolant which is not too far from the optimal polynomial P Opt N , i.e.

f -P Tch N ∞ 1 + Λ Tch N f -P Opt N ∞ , * Carl Gustav Jacob Jacobi (1804 -1851
) is a German Mathematician. † Jacques Charles François Sturm (1803 -1855) is a French Mathematician who was born in Geneva which was a part of France at that time.

‡ Joseph Liouville (1809 -1882) is a French Mathematician who founded the Journal de Math ḿatiques Pures et Appliquées.

§ Adrien-Marie Legendre (1752 -1833) is a French Mathematician. Not to be confused with an obscure politician Louis Legendre. Their portraits are often confused.

¶ Pafnuty Tchebyshev (1821 -1894), a Russian mathematician who contributed to many fields of Mathematics from number theory to probabilities and numerical analysis. where f : [-1, 1] → R is the function that we interpolate and P Tch N is the interpolation polynomial constructed on nodes (2.5). Here Λ Tch N is the so-called Lebesgue * constant for Tchebyshev interpolation polynomials of degree N. Notice that Λ Tch N does not depend on the function f being interpolated. The Lebesgue constants for various interpolation techniques have the following asymptotic behaviour

Λ Tch N ∼ O(log N) , Λ Leg N ∼ O( √ N) , Λ Uni N ∼ O 2 N N log N ,
where Λ Leg N and Λ Uni N are Lebesgue constants for Legendre and uniform nodes distributions correspondingly. In particular one can see that the uniform node distribution is simply disastrous. This phenomenon is another 'avatar' of the so-called Runge phenomenon illustrated in Figure 2. The performance of Tchebyshev's nodes is shown in Figure 3. * Henri Lebesgue (1875 -1941), a French mathematician most known for the theory of integration having its name. Remark 3. As it was shown by [START_REF] Vértesi | Optimal Lebesgue constant for Lagrange interpolation[END_REF] [39], the Lebesgue constant for Tchebyshev distribution of nodes is very close for the smallest possible Lebesgue constant:

Λ Tch N = 2 π ln N + γ + ln 8 π + o (1), Λ min N = 2 π ln N + γ + ln 4 π + o (1) .
where γ ≈ 0.5772156649 . . . is the Euler * -Mascheroni † constant.

The Lebesgue constant. It is worth to explain better the important notion of the Lebesgue constant and how it appears in the theory of interpolation (below we follow [START_REF] Solin | Partial Differential Equations and the Finite Element Method[END_REF]). Let us take an arbitrary (valid) nodes distribution {x i } n i = 0 in domain U, i.e. ∀i = 1, 2, . . . , n :

x i ∈ U , ∀j = i : x i = x j .
For any continuous function u ∈ C(U) there exists a unique interpolating polynomial P(x) ∈ P n [ R ] of degree n = deg P :

P n (x i ) = u i ≡ u (x i ) , i = 1, 2, . . . , n .
This interpolating polynomial P n (x) will be denoted as I n [ u ] . We introduce also the best possible approximating polynomial P * ∈ P n [ R ] :

u -P * ∞ = inf P ∈ Pn [ R ] u -P ∞ .
It is not obliged that P * ≡ I n [ u ] . However, since P * ∈ P n [ R ] we necessarily have

P * ≡ I n [ u ]
. Therefore, the following inequalities hold:

u -I n [ u ] ∞ = u -P * + P * -I n [ u ] ∞ ≡ u -P * + I n [ P * ] -I n [ u ] ∞ u -P * ∞ + I n [ P * ] -I n [ u ] ∞ u -P * ∞ + I n • P * -u ∞ 1 + I n • u -P * ∞ .
The norm of a linear operator I n [•] is defined as

I n def := sup u ∞ = 1 I n [ u ] ∞ .
The aforementioned norm of the interpolation operator I n [•] is called the Lebesgue constant for the set of nodes {x i } n i = 0 . In multiple dimensions the Lebesgue constant depends also on the shape of domain U additionally to the nodes distribution {x i } n i = 0 . As we mentioned above, obvious choices such as the uniform distribution of nodes is disastrous, since it yields the exponential growth of the Lebesgue constant Λ Uni n ∼ O 2 n n log n , where n is the degree of the interpolating polynomial.

The estimation of the Lebesgue constant in multi-dimensional non-Cartesian domains is a problem essentially open nowadays. The most precious information is to find the nodes distribution, for example over a triangle, which minimizes the Lebesgue constant. This knowledge would be crucial for the design of new spectral elements [START_REF] Solin | Partial Differential Equations and the Finite Element Method[END_REF]. The locations of nodes which minimize the magnitude of the Lebesgue constant Λ n → min for a given polynomial space P n [ R ] are called Lebesgue points. Nowadays, the best known points on quadrangles are tensor products of Gauß-Lobatto * or Tchebyshev nodes. On triangles the Fekete † points seem currently to be the best choice.

Intermediate conclusions. The remarks above show that Tchebyshev polynomials satisfy the requirement (1) for a useful spectral basis. The requirement (2) is met due to derivative recursion formulas, which can be easily demonstrated (see also Appendix A):

T ′ k (x) = 2 k [ k-1 2 ] n = 0 1 δ k-1-2 n T k-1-2 n (x) , δ k = 2 , k = 0 , 1 , k = 0 .
Finally, the requirement (3) is satisfied as well thanks to the Fast Cosine Fourier Transform (FCFT) (a variant of FFT). It allows to compute spectral coefficients {v k } N k = 0 from the node values u n (x i ) N i = 0 and vice versa. Consequently, Tchebyshev polynomials have become an almost universal choice for non-periodic problems. These methods in 1D have been implemented in the Matlab toolbox Chebfun * (and Chebfun2 in 2D).

Remark 4. To Author's knowledge, Legendre polynomials found some applications in the construction of Spectral Element Method (SEM) bases [START_REF] Solin | Partial Differential Equations and the Finite Element Method[END_REF].

Remark 5. We have to mention that high-order polynomial interpolants have a bad reputation in the Numerical Analysis (NA) community for the two following reasons:

(1) Because of the Runge phenomenon (2) and because of the following Theorem 1. For any node density distribution function there exists a continuous function such that the L ∞ -norm of the interpolation error tends to infinity when the number of nodes N → +∞.

However, this bad reputation is not justified. For instance, the Runge phenomenon is completely suppressed by Tchebyshev nodes distribution. The low or moderate Lebesgue constant Λ N ensures that we are not too far from the optimal polynomial.

Infinite domains. The infinite domains can be handled, first of all, by periodisation in the context of Fourier-type methods as explained above. However, truly infinite domains can be handled by various approaches:

• Hermite polynomials • Sinc functions • (almost) Rational functions • Change of variables, e.g. x : [-π, π] → R , x(q) = ℓ tan q 2 , ℓ ∈ R + • Truncation of a domain (-∞ , ∞) [-ℓ, ℓ]
followed by a change of variables such as

x : [-π, π] → [-ℓ, ℓ] ,
x(q) = ℓ q π . * Chebfun's team is lead by Nick Trefethen.

Two last items are followed by a conventional pseudo-spectral discretization on a finite domain [-π, π].

A Hermite * -type pseudo-spectral method can be briefly sketched as follows. First of all, we construct recursively the family of Hermite's polynomials:

H 0 (x) = 1 , H 1 (x) = 2 x, H n+1 (x) = 2 x H n (x) -2 n H n-1 (x), n 1 .
Then, we construct Hermite's functions, which form an orthonormal basis on R:

H n (x) def := 1 √ 2 n • n! H n (x) e -x 2 2 .
Derivatives of expansions in Hermite functions can be easily re-expanded using the relation:

H ′ n (x) = - n + 1 2 H n+1 (x) + n 2 H n-1 (x) .
However, the lack of a Fast Hermite Transform (FHT) implies the use of differentiation matrices in numerical implementations. There are also some concerns about a poor convergence rate when the number of modes N is increased. The infinite domains remain very challenging in silico.

Determining expansion coefficients

There are three main techniques to find the spectral expansion coefficients {v k } N k=0 . In order to explain them, let us introduce first the residual function on a trial solution u n (x, t):

R[u n ](x, t) def := L u n -g (x, t).
Typically, we shall evaluate the residual R on the expansion (2.2) and the residual norm R is generally considered as a measure of the approximate solution quality. The goal is to keep the residual as small as possible across the domain U (in our particular case

U = [-1, 1]).
So, we can mention here at least five approaches to determine the expansion coefficients:

Tau-Lanczos: Spectral coefficients {v k } N k=0 are selected such that the boundary conditions are satisfied identically and the residual R[u n ] is orthogonal to as many basis functions φ k (x) as possible.

Galerkin: First the basis functions are recombined φ k (x) N k=0 φk (x) N k=0 so that the boundary conditions are satisfied identically. Then, the coefficients {v k } N k=0 are * Charles Hermite (1822 -1901), a French mathematician whose "Cours d'analyse" represent a lot of interest even today. He was also the PhD advisor of another great French mathematician -Henri Poincaré.

found so that the residual R[u n ] be orthogonal to as many of new basis functions φk (x) N k=0 as possible.

Collocation: This approach is similar to the Tau-Lanczos method concerning the boundary conditions: spectral coefficients {v k } N k=0 are selected such that the boundary conditions are satisfied. The rest of coefficients is determined so that the residual R[u n ](x, t) vanishes at as many (thoroughly chosen) spatial locations as possible.

Petrov-Galerkin: It is a variant of Galerkin method in which the residual R[u n ] is made orthogonal to a set of functions, which is different from the approximation space basis φ k (x) N k=0 .

Least squares: Various least square-type approaches are used when for some reason the number of coefficients to be determined is different from the number of conditions which can be imposed. Below we shall avoid such pathological situations.

The Tau-Lanczos technique was proposed by C. Lanczos * in 1938. What we call the Galerkin † technique was proposed independently first by I. Bubnov ‡ and by W. Ritz § (one more example of the Arnold principle in action!). So, to respect the historical time line, the method should be called Bubnov-Ritz-Galerkin. Today it is the basis of the Finite Element Method (FEM). The collocation technique was called the pseudo-spectral method presumably for the first time by S. Orszag ¶ in 1972. The Petrov-Galerkin method was proposed by G.I. Petrov . It is used up to now in some convection-dominated problems.

Remark 6. The collocation method can be recast into the Petrov-Galerkin framework when we make the residual R[u n ] is made orthogonal to Dirac * * singular measures δ(xx k ) n k=1 , where {x k } n k=1 are the collocation points. * Cornelius Lanczos (1893 -1974), a Hungarian/American numerical analyst. His books are very pedagogical as well.

† Boris Galerkin (more precisely romanized as Galyorkin) (1871 -1945), a Russian then Soviet Civil Engineer. The Author suggests to read his biography which is comparable to James Bond (007) movies.

‡ Ivan Bubnov (1872 -1919), a Russian naval engineer.

§ Walther Ritz (1878 -1909), a talented Swiss Physicist, who died young in Göttingen, Germany.

¶ Steven Orszag (1943 -2011), an American numerical analyst, one of the first users of pseudo-spectral methods.

Georgii Ivanovich Petrov (1912 -1987), a Russian fluid mechanician. * * Paul Adrien Maurice Dirac (1902 -1984) is a British Theoretical Physicist who predicted theoretically using Dirac's equation the existence of the positron. One of the founders of the Quantum Mechanics.

For linear problems all these methods work equally well. However, for nonlinear ones (and in the presence of variable coefficients) the pseudo-spectral (collocation) approach is particularly easy to apply since it involves the products of numbers (solution/variable coefficient's values in collocation points) instead of products of expansions, which are much more difficult to handle.

The convergence of pseudo-spectral approximations for very smooth functions is always geometrical, i.e. ∼ O(q N ), where N is the number of modes. This statement is true for any derivative with the same convergence factor 0 < q < 1. However, the periodic pseudo-spectral method converges always faster than its non-periodic counterpart. This conclusion follows from convergence properties of Fourier and Tchebyshev series in the complex domain.

The relative resolution ability of various pseudo-spectral methods can be also quantified in terms of the number of points per wavelength needed to resolve a signal. Indeed, this description is more suitable for wave propagation problems. For periodic Fourier-type methods one needs 2 points per wavelength. For Tchebyshev-type methods one needs about π points. Finally, this number goes up to 6 nodes per wavelength for uniform grids. However, due to the huge Lebesgue constant Λ Uni N the uniform grids in pseudo-spectral setting are not usable in practice.

Aliasing, interpolation and truncation

Let us take a continuous (and possibly a smooth) function u(x) defined on the interval I = (-π, π) and develop it in a Fourier series. In general it will contain the whole spectrum (i.e. the infinite number) of frequencies:

u(x) = +∞ k = -∞ v k e ikx .
Now let us discretize the interval I with N equispaced collocation points (as we do it in Fourier-type pseudo-spectral methods. In the following we assume N to be odd, i.e. N = 2m + 1. On this discrete grid all modes e i(k+jN )x j∈Z are indistinguishable. See Figure 4 for an illustration of this phenomenon.

The interpolating trigonometric polynomial on a given grid can be written as

I N [u] = m k = -m vk e ikx .
Each discrete Fourier coefficient incorporates the contributions of all modes which looks the same on the considered grid: Let us recall that the polynomial

vk = +∞ j = -∞ v k+jN . x -3 -2 -1 0 1 2 3 y(x)
I N [u] takes the prescribed values u(x k ) m k = -m in the points of the grid {x k } m k = -m .
This object is fundamentally different from the truncated Fourier series:

T N [u] = m k = -m v k e ikx .
The difference between these two quantities is known as the aliasing error :

R N [u] def := I N [u] -T N [u] = m k = -m ∀j =0 v k+jN e ikx .
After applying the Pythagoras theorem * , we obtain

u -I N [u] 2 L 2 = u -T N [u] 2 L 2 + R N [u] 2 L 2 .
Thus, the interpolation error is always larger than the truncation error in the standard L 2 norm. The amount of this difference is precisely equal to the committed aliasing error. However, we prefer to use in pseudo-spectral methods the interpolation technique because of the Discrete Fourier Transform (DFT), which allows to transform quickly (thanks to the FFT algorithm) from the set of function values in grid points to the set of its interpolation coefficients. So, it is easier to apply an FFT instead of computing N integrals to determine the Fourier series coefficients.

Nonlinearities. Let us take the simplest possible nonlinearity -the product of two functions u(x) and v(x) defined by their truncated Fourier series containing the modes up to m :

u(x) = m k = -m u k e ikx , v(x) = m k = -m v k e ikx .
The product of these two functions w(x) can be obtained my multiplying the Fourier series:

w(x) = u(x) • v(x) = m k = -m u k e ikx • m k = -m v k e ikx ≡ 2 m k = -2 m w k e ikx .
It can be clearly seen that the product contains high order harmonics up to e ±imx which cannot be represented on the initial grid. Thus, they will contribute to the aliasing error explained above. The aliasing of a nonlinear product can be ingeniously avoided by adopting the so-called 3/2 th rule whose Matlab implementation is given below. This function assumes that input vectors are Fourier coefficients of functions u(x) and v(x) . The resulting vector contains (anti-aliased) Fourier coefficients of their product w(x) = u(x) • v(x) . The main idea behind is to complete vectors of Fourier coefficients by a sufficient number of zeros (i.e. the so-called zero padding technique) so that in the physical space the product u(x) • v(x) can be fully resolved. The final step consists in extracting m relevant Fourier coefficients [START_REF] Trefethen | Spectral methods in MatLab[END_REF][START_REF] Uecker | A short ad hoc introduction to spectral methods for parabolic PDE and the Navier-Stokes equations[END_REF].

Remark 7. To Author's knowledge, the development of efficient and rigorously justified anti-aliasing rules for other types of nonlinearities such as the division, square root, etc. is an open problem.

Example of a second order boundary value problem

Consider the following second order Boundary Value Problem (BVP) on the interval I = [-1, 1]:

u xx + u x -2u + 2 = 0 , u(-1) = u(1) = 0 . (3.1)
This BVP has the exact solution

u(x) = 1 - sinh(2) sinh(3) e x - sinh(1) sinh(3) e -2x , x ∈ [-1, 1] . (3.2)
We shall seek for the numerical solution to (3.1) in the form of a truncated Tchebyshev expansion:

u(x) ≈ 4 k = 0 a k T k (x) , ∀k : a k ∈ R Spectral coefficients {a k } 4 k=0
have to be determined. The last expansion is substituted into the governing equation (3.1). There is no reason that (3.1) will be satisfied identically in every point of I. Hence, we can measure the residual

R(x) = u xx + u x -2u + 2 (x) 0 .
The enforcing of boundary conditions u(-1) = u(1) = 0 leads to two additional relations on spectral coefficients:

a 0 + a 1 + a 2 + a 3 + a 4 = 0 , a 0 -a 1 + a 2 -a 3 + a 4 = 0 .
So, we have two relations coming from boundary conditions and we have five degrees of freedom {a k } 4 k=0 . It means that we have to impose three additional conditions to determine uniquely all spectral coefficients. Different approaches prescribe different numerical recipes.

Tau-Lanczos

We require that the residual R(x) be orthogonal to the first three basis functions T 0,1,2 (x). It gives us three additional relations:

R, T k ≡ ˆ1 -1 R(x) T k (x) √ 1 -x 2 dx = 0 , k = 0, 1, 2 .

Galerkin

We recombine the Tchebyshev polynomials to form a different basis which satisfies identically the boundary conditions: Then we require that the residual R(x) is orthogonal to the new basis functions:

φ 0 (x) def := T 2 -T 0 (x) , φ 1 (x) def := T 3 -T 1 (x) , φ 2 (x) def := T 4 -T 0 (x) .
R, φ k ≡ ˆ1 -1 R(x) φ k (x) √ 1 -x 2 dx = 0 , k = 0, 1, 2 .

Collocation

This approach is particularly simple. We require that R(x k ) ≡ 0 in three interior Tchebyshev points x k = cos πk 4 , k = 1, 2, 3.

Conclusions. The comparison of three numerical solutions with the exact one (3.2) is shown in Figure 5. In this Figure 5 one can see that the collocation * and Galerkin methods provide nearly identical † numerical solutions which are closer to the exact solution than the prediction of the Tau method. However, in general the Author is impressed by the performance of spectral methods -with only five degrees of freedom {a k } 4 k=0 we capture very well the global behaviour of the solution (3.2) in every point of the interval I. Perhaps, the last point has to be emphasized again: the (pseudo-)spectral methods provide the numerical value of the solution in the whole domain, not only in collocation or grid points.

Application to heat conduction

In this Section we show how to apply the Fourier-type spectral methods to simulate some simple and not so simple diffusion processes.

An elementary example

As the simplest example, consider the linear heat equation

u t = ν u xx , x ∈ R , (4.1)
completed with the initial condition u(x, 0) = u 0 (x) .

We consider this equation on the whole line R for the sake of simplicity. However, in silico it will be truncated to a periodic interval (in the example below it will be

I = [-1, 1] ).
Let us apply the Fourier integral transform to the both sides of equation 4.1:

dû(k, t) dt = -ν k 2 û(k, t) , (4.2)
where the forward û(k, t) = F{u(x, t)} and inverse u(x, t) = F -1 {û(k, t)} integral Fourier transforms are defined below in (4.7) and (4.8) correspondingly. The transformed heat equation (4.2) can be regarded as a linear ODE. Its solution can be readily obtained:

û(k, t) = û0 (k) e -ν k 2 t , û0 (k) = F{u 0 (x)} .
Consequently, we possess an analytical solution to the heat equation (4.1) in Fourier space. In order to obtain the solution in physical space, an inverse Fourier transform has to be computed. The exponential decay of Fourier coefficients û(k, t) (except the zero mode k = 0) ensures that the solution u(x, t) becomes infinitely smooth C ∞ for t > 0 . For example, the initial condition and corresponding solution at t = T = 5 s are represented in Figure 6. The Matlab code used to generate Figure 6 is provided below as well.

A less elementary example

The numerical code above is based on the knowledge of the analytical solution to ODE (4.2) in the Fourier space. It is unnecessary to say that an analytical solution is available in very simple situations only. Consequently, we assume that the resulting ODE system in the physical (Fourier) space is formally solved (i.e. advanced in time) by the semigroup operator

S(t) ( Ŝ(t) ≡ F • S(t) • F -1 correspondingly): u(x, t) ≡ S(t) • u 0 (x) .
In practice this semigroup operator is realized using numerical time-marching techniques described briefly in Appendix F. So, a more general Fourier spectral algorithm is (1) Decompose the initial condition:

û0 (k) = F u 0 (x) , (2) Advance in time: ûk (t) = Ŝ(t) û0 (k) , (3) Synthesize: u(x, t) = F -1 ûk (t)
. This algorithm is based on the fact that the following diagram commutes:

u(x, 0) u(x, t) û0 (k) ûk (t) F S(t) Ŝ(t)
F -1

A real-life example

In this Section we shall consider a realistic model (to be honest we take a slightly simplified version), which was proposed in [START_REF] Mendes | A method for predicting heat and moisture transfer through multilayered walls based on temperature and moisture content gradients[END_REF] to predict heat and moisture transfer through the walls. This model is used in real-world Civil Engineering applications such as the code Domus. So, the model we consider in this Section reads (for simplicity we consider the 1D case):

∂θ ∂t = ∂ ∂x D θ ∂θ ∂x + D T ∂T ∂x , (4.3
)

ρc m ∂T ∂t = - ∂q ∂x -L(T ) • ∂j v ∂x , (4.4)
where θ is moisture volumetric content (i.e. moisture density) and T is the temperature. The mass transport coefficients D θ (θ, T ) and D T (θ, T ) may depend nonlinearly on the solution θ(x, t), T (x, t) . Here we assume for simplicity that the mass density ρ and the specific moisture heat c m are some positive constants * . Finally, the heat flux q and vapor flow j v can be expressed as

q = -λ(θ, T ) ∂T ∂x , j v = -V θ (θ, T ) ∂θ ∂x -V T (θ, T ) ∂T ∂x .
In Fourier-type pseudo-spectral methods we usually work with Fourier coefficients. Consequently, we apply the Fourier transform to both sides of equations (4.3), (4.4):

d θ dt = ik F D θ ∂θ ∂x + D T ∂T ∂x , (4.5
)

ρc m d T dt = ik F λ(θ, T ) ∂T ∂x -F L(T ) • ∂j v ∂x , (4.6)
where we introduced some notations for the Fourier transform F(•):

θ(k, t) def := F θ(x, t) = ˆ+∞ -∞ θ(x, t) e ikx dx , (4.7)
where k ∈ R is the wave number. Inversely we have

θ(x, t) = F -1 θ(k, t) = 1 2π ˆ+∞ -∞ θ(k, t) e -ikx dk . (4.8)
Now, system (4.5), (4.6) can be considered as a coupled nonlinear system of Ordinary Differential Equations (ODEs) (not PDEs!) for the Fourier coefficients of solutions θ(x, t), T (x, t) . The numerical methods for systems of ODEs will be explained in a separate course. As general fundamental references on this topic we can recommend [START_REF] Hairer | Geometric Numerical Integration[END_REF][START_REF] Hairer | Solving ordinary differential equations: Nonstiff problems[END_REF]. The Author even suggests to employ a well-documented ready-to-use ODE library such as [START_REF] Shampine | The MATLAB ODE Suite[END_REF], for example. Now we have to explain how to evaluate the right hand side of equations (4.5), (4.6) when only the Fourier coefficients are available. The recipe is very simple:

• All linear operations (e.g. additions, subtractions, multiplication by a scalar) can be equally made in Fourier or in a real spaces. The choice has to be done in order to minimize the number of FFT operations

F αθ + βT ≡ α θ + β T
• All nonlinear products are made in the real space and then we transfer the result back to the Fourier space using one FFT, e.g.

F θ(x, t) • T (x, t) = F F -1 θ(k, t) • F -1 T (k, t)
• All spatial derivatives are computed only in Fourier space, as

∂ n ∂x n • ⇆ F -1 (ik) n •
The Reader can notice that above we used Fourier integrals. It comes from the mathematical tradition. In computer implementations one has to take a finite interval, periodize it and use direct and inverse Discrete Fourier Transforms (DFTs) instead of F{•} and

F -1 {•} respectively.
The computation of derivatives using the Fourier collocation spectral method is illustrated below on the following periodic function (with period 2):

u(x) = sin π(x + 1) e sin(π(x + 1)) , x ∈ [-1, 1] . (4.9)
The first three derivatives of this function can be readily computed using any symbolic computation software (the Author used Maple):

u ′ (x) = π cos π(x + 1) 1 + sin π(x + 1) e sin(π(x + 1)) , u ′′ (x) = π 2 cos 2 π(x + 1) sin π(x + 1) + 3 -sin π(x + 1) -1 e sin(π(x + 1)) , u ′′′ (x) = π 3 cos π(x + 1) cos 2 π(x + 1) sin π(x + 1) + 6 -7 sin π(x + 1) -4 .

The expressions above are used to assess the accuracy of the computed approximations in collocation points. The analytical derivatives with computed ones (red dots) are shown in Figure 7 for N = 128 collocation points. To the graphical accuracy the results are indistinguishable. That is why we computed also the discrete ℓ ∞ norms to compute the relative errors:

ε (p) N = u num -u exact ∞ u exact (x i ) ∞ .
For N = 32 Fourier modes we obtain the following numerical results:

ε (1) 32 ≈ 1.5 × 10 -15 , ε (2) 32 ≈ 8.4 × 10 -15 , ε (3) 32 ≈ 4.7 × 10 -14 .
The Matlab code used to generate these results and Figure 7 is provided below. Above we exploit the property that the Fourier transform is a 'change a variables' where the differentiation operator ∂ x becomes diagonal. Fourier collocation spectral method. We use 128 collocation points. The agreement is perfect up to graphic accuracy.

Indications for further reading

If you got interested in the beautiful topic of pseudo-spectral methods, you can find more information in the following books: • The book by J. Boyd * is probably the most exhaustive one. It covers many topics and applications of spectral methods. The material is presented as a collection of tricks. For example, it is one of seldom books where (semi-)infinite domains and spherical geometries are covered. I am not sure that after reading this book you will know how to program the pseudo-spectral methods, but you will have a broad of view of possible issues and how to address them [

A. Some identities involving Tchebyshev polynomials

Tchebyshev polynomials {T n (x)} ∞ n=0 form a weighted orthogonal system on the segment [-1, 1]:

T n , T m ≡ ˆ1 -1 T n (x) • T m (x) √ 1 -x 2 dx =        0, m = n , π, m = n = 0 , π 2 , m = n > 0 .
The first few Tchebyshev polynomials are

T 0 (x) = 1 , T 1 (x) = x , T 2 (x) = 2 x 2 -1 , T 3 (x) = 4 x 3 -3 x , T 4 (x) = 8 x 4 -8 x 2 + 1 , T 5 (x) = 16 x 5 -20 x 3 + 5 x .
They are represented graphically in Figure 8. Higher order Tchebyshev polynomials can be constructed using the three-term recursion relation:

T n+1 (x) = 2 x T n (x) -T n-1 (x) .
The following Matlab code realizes this idea in practice:

1 function P = Chebyshev(n) Then, n th Tchebyshev polynomial can be evaluated using the standard Matlab's function polyval(), e.g. There is an explicit expression for the n th polynomial:

T n (x) = cos nθ(x) , θ(x) def := arccos x . (A.1)
The first derivatives of Tchebyshev polynomials can be constructed recursively as well:

T ′ n+1 (x) n + 1 = T ′ n-1 (x) n -1 + 2 T n (x) . (A.2)
Otherwise, the first derivative T ′ n (x) can be found from the following relation:

(1 -x 2 ) T ′ n (x) = -n x T n (x) + n T n-1 (x)
. Tchebyshev polynomial T n (x) satisfies the following linear second order differential equation with non-constant coefficients:

(1 -x 2 ) T ′′ n -x T ′ n + n 2 T n = 0 .
Finally, zeros of the n th Tchebyshev polynomial T n (x) are located at

x 0 k = cos 2k -1 2n π , k = 1, 2, . . . , n,
and extrema at

x ext k = cos πk n , k = 0, 1, . . . , n .
When implementing Tchebyshev-type pseudo-spectral methods for nonlinear problems, one can use the expression of the product of two Tchebyshev polynomials in terms of higher and lower degree polynomials:

T m (x) T n (x) = 1 2 T n+m (x) + T n-m (x) , ∀n m 0 .
To finish this Appendix we give the generating function for Tchebyshev polynomials:

1 -x z 1 -2 x z + z 2 = ∞ n = 0 T n (x) z n .
In order to satisfy the requirement (2) from Section 2.1, we provide here the relations which allow to re-express the derivatives of the Tchebyshev expansion in terms of Tchebyshev polynomials again. Namely, consider a truncated series expansion of a function u(x) in Tchebyshev polynomials:

u(x) = N k = 0 v k T k (x) .
Imagine that we want to compute the first derivative of the expansion above (it is always possible, since the sum is finite and we can differentiate term by term):

u ′ (x) = N k = 0 v k dT k (x) dx .
Now we re-expand the derivative u ′ (x) in the same basis functions:

u ′ (x) = N k = 0 v ′ k T k (x) .
The main question is how to re-express the coefficients {v ′ k } N k=0 in terms of coefficients {v k } N k=0 ? This goal is achieved using basically the recurrence relation (A.2) (even if it does not jump into the eyes). Thanks to it we have an explicit relation between the coefficients:

v ′ k = 2 δ k N j = k+1 j+k odd j v j , j = 0, 1, . . . , N -1 ,
where

δ k = 2, k = 0 , 0, k = 0 .
A similar 'trick' can be made for the 2 nd derivative as well:

u ′′ (x) = N k=0 v k d 2 T k (x) dx 2 = N k=0 v ′′ k T k (x) .
The connection between coefficients {v ′′ k } N k=0 and {v k } N k=0 is given by the following explicit formula:

v ′′ k = 1 δ k N j=k+2 j+k even j (j 2 -k 2 ) v j , j = 0, 1, . . . , N -2 .
Finally, in order to construct the spectral coefficients for the n th derivative, one can use the following recurrence relation (also stemming from (A.2)):

δ k-1 v (n) k-1 = v (n) k+1 + 2 k v (n-1) k , k 1 ,
which has to be completed by starting value v ′ N ≡ 0 . Remark 8. Similar identities exist for the family of Jacobi polynomials P α, β n (x) as well. However, they are more complicated due to the presence of two arbitrary parameters α, β > -1. The family of Tchebyshev polynomials is a particular case of Jacobi polynomials when α = β = -1 2 .

A.1. Compositions of Tchebyshev polynomials

We would like to show here another interesting formula involving Tchebyshev polynomials:

Theorem 2 (Composition formula). If T n (x) and T m (x) are Tchebyshev polynomials (m, n 0), then T m • T n (x) ≡ T m T n (x) = T mn (x) . (A.3)
Proof. An eleven (!) pages combinatorial proof of this result can be found in [START_REF] Benjamin | Combinatorially composing Chebyshev polynomials[END_REF]. Here we shall prove the composition formula (A.3) in one line by using some rudiments of complex variables. Let us introduce the variable z def := e iθ ∈ * such that

x = 1 2 z + 1 z ⇆ x = cos θ .
Then, by (A.1) we have a complex representation of the n th Tchebyshev polynomial:

T n (x) ≡ T n 1 2 z + z -1 = 1 2 z n + 1 z n
. Now, we have all the elements to show the main result:

T m T n 1 2 z + z -1 = T m 1 2 z n + z -n = 1 2 (z n ) m + 1 (z n ) m = 1 2 z mn + 1 z mn = T mn 1 2 z + z -1 ≡ T mn (x) .

B. Trefftz method

In this Appendix I would like to give the flavour of the so-called Trefftz † methods, which remain essentially unknown/forgotten nowadays. These methods belong to * Symbol denotes the unit circle S 1 (0) on the complex plane C. † Erich Emmanuel Trefftz (1888Trefftz ( -1937)), a German Mathematician and Mechanical Engineer.

boundary-type solution procedures. Below I shall quote a Professor * from Laboratoire Jacques-Louis Lions (LJLL) at Paris 6 Pierre and Marie Curie University, who makes interesting comments about the knowledge of Trefftz methods in French Applied Mathematics community (the original language and orthography are conserved):

[ . . . ] Sinon j'ai moi aussi fait une petite enquête, à la suite de laquelle il apparaît que personne ne connait Trefftz chez les matheux de Paris 6, hormis Nataf qui en a entendu parler pendant son DEA en méca !!! Trefftz est un grand oublié car son papier est vraiment extrêmement intéressant, encore plus si tu te rends compte qu'il produit en fait une estimation a posteriori (ça doit même être la première). En revanche je me suis aussi persuadé que les méthodes de Trefftz végètent chez les mécano. [ . . . ] So, the situation is rather sad. Fortunately this method continued to live in the Mechanical Engineering community who preserved it for future generations.

In the exposition below we shall follow an excellent review paper [START_REF] Kita | Trefftz method: an overview[END_REF]. More precisely we describe the so-called indirect Trefftz method, which was proposed originally by E. [START_REF] Trefftz | Gegenstück zum ritzschen Verfahren[END_REF] in [START_REF] Trefftz | Gegenstück zum ritzschen Verfahren[END_REF]. There are also direct Trefftz methods proposed some sixty years later in [4]. They are much closer to the Boundary Integral Equation Methods (BIEM) and will not be covered here. Interested readers can refer to [START_REF] Kita | Trefftz method: an overview[END_REF].

Consider a (compact) domain Ω ⊆ R d , d 2 and a (scalar) equation on it (see Figure 9 for an illustration):

L u = 0 . (B.1)
You can think, for example, that

L = -∇ 2 = -d i = 1 ∂ 2 ∂x 2 i
is the Laplace † operator.

Equation (B.1) has to be completed by appropriate boundary conditions:

B u = u • , x ∈ ∂ Ω , (B.2)
where u • is the boundary data (solution value or its flux) and B is an operator depending on the type of boundary conditions in use. For example, in the case of Laplace equation the following choices are popular:

Dirichlet: B def := I (Identity operator, i.e. Iu ≡ u); in this case the solution value is prescribed on the boundary.

Neumann:

B def := ∂ ∂n ≡ n • ∇ = d i = 1 n i • ∂ ∂x i
(Normal derivative, n being the exterior normal to ∂ Ω); physically it corresponds to the prescribed flux through the boundary. * His/Her identity will be hidden in order to avoid any kind of diplomatic incidents. † Pierre-Simon de Laplace (1749 -1827) is a French Mathematician whose works were greatly disregarded during his times. The understanding of their importance came much later (il vaut mieux tard que jamais ).

Robin: B

def := α I + β ∂ ∂n
, where α, β ∈ R are some parameters; this case is a mixture of two previous situations. It arises in some problems as well (see e.g. [START_REF] Chhay | Asymptotic heat transfer model in thin liquid films[END_REF] for the heat conduction problem in thin liquid films).

Notice, that the boundary ∂ Ω can be divided in some sub-domains where a different boundary condition is imposed:

∂ Ω = Γ 1 ∪ Γ 2 ∪ . . . ∪ Γ n µ d Γ i ∩ Γ j = 0 , 1 i < j n ,
where µ d is the Lebesgue measure in R d .

In the indirect Trefftz method the numerical solution is sought as a linear combination of T-complete functions φ k (x) [START_REF] Herrera | Boundary Methods: An Algebraic Theory[END_REF], which satisfy exactly the governing equation (B.1):

n k = 1
u n (x) = n k = 1 v k φ k (x) . (B.3)
For example, for the Laplace equation in R 2 the T-complete set of functions is [START_REF] Herrera | Boundary methods: development of complete systems of solutions[END_REF] r k e ikθ ∞ k=0 , where (r, θ) are polar coordinates * on plane. The coefficients {v k } n k=1 are to be determined. To achieve this goal the approximate solution (B.3) is substituted into the boundary conditions (B.2) to form the residual:

R[u n ] = B u n -u • .
If the residual R[u n ] is equal identically to zero on the boundary ∂ Ω, then we found the exact solution (we are lucky!). Otherwise (we are unlucky and in Numerical Analysis it happens more often), the residual has to be minimized. Very often the boundary operator B is linear and we can write:

R[u n ] = n k = 1 v k B φ k (x) -u • .
So, we can apply now the collocation, Galerkin or least square methods to determine the coefficients {v k } n k = 1 . These procedures were explained above. Remark 9. Suppose that a T-complete set of functions for equation (B.1) is unknown. However, there is a way to overcome this difficulty if we know the Green † function ‡ * r =

x 2 + y 2 , θ = arctan y x , with arctan (±∞) = ± π 2 . † George Green (1793 -1841), a British Mathematician who made great contributions to the Mathematical Physics and PDEs.

‡ The Green function is named after George Green (see the footnote above), but the notion of this function can be already found first in the works of Laplace, then in the works of Poisson. As I said, the works of Laplace were essentially disregarded by his "colleagues". Poisson was luckier in this respect. He is also known for (mis)using his administrative resource in order to delay (just for a couple of decades, nothing serious ) the publication of his competitors, e.g. the young (at that time) Cauchy. G(x; Q), x ∈ Ω for our equation (B.1), i.e.

L G = δ x -Q , x ∈ Ω , Q ∈ R d .
where δ(x) is the singular Dirac measure. Notice that point Q can be inside or outside of domain Ω, where the problem is defined. Then, we can seek for an approximate solution u n (x) as a linear combination of Green functions:

u n (x) = n k = 1 v k G k (x; Q k ) ,
where the points Q k n k=1 are distributed outside of the computational domain Ω in order to avoid the singularities. This method is schematically illustrated in Figure 9. The accuracy of this approach depends naturally on the distribution of points Q k n k=1 . To Author's knowledge there exist no theoretical indications for the optimal distribution (as it is the case of Tchebyshev nodes on a segment). So, it remains mainly an experimental area of the research.

Conclusions. In the indirect Trefftz methods (historically, the original one) the problem solution is sought as a linear combination of the functions, which satisfy the governing equation identically. Then, the unknown coefficients are chosen so that the approximate solution satisfies the boundary condition(s) as well. It can be done by means of collocation, Galerkin or least square procedures. The main advantage of Trefftz methods is that it allows to reduce problem's dimension by one (i.e. 3D

2D and 2D 1D), since the residual R[u n ] is minimized on the domain boundary ∂ Ω ⊆ R d-1 . Readers who are interested in the application of Trefftz methods to their problems should refer to [START_REF] Kita | Trefftz method: an overview[END_REF] and references therein.

C. A brief history of diffusion in Physics

Since the main focus of the PhD school is set on the diffusion processes (molecular diffusion, heat and moisture conduction through the walls, etc. ), it is desirable to explain how this research started and why the diffusion is generally modeled by parabolic PDEs [START_REF] Evans | Partial Differential Equations[END_REF]. The historic part of this Appendix is partially based on [START_REF] Philibert | One and a half century of diffusion: Fick, Einstein, before and beyond[END_REF].

Firs of all, let us make a general, perhaps surprising, remark, which is based on the analysis of historical investigations: understanding the microscopic world is not compulsory to propose a reliable macroscopic law. Indeed, Fourier did not know anything about the nature of heat, Ohm * about the nature of electricity, Fick † about salt solutions and Darcy ‡ about the structure of porosity and water therein. In cases of the diffusion, the * Georg Simon Ohm (1789 -1854) was a German Physicist.

† Adolf Eugen Fick (1829 -1901) was a German physician and physiologist. ‡ Henry Darcy (1803 -1858) was a French engineer in Hydraulics. bridge between microscopic and macroscopic worlds was built by A. Einstein * in 1905. Namely, he expressed a macroscopic quantity -the diffusion coefficient -in terms of microscopic data (this result will be given below). In general, 1905 was Annus Mirabilis † for A. Einstein. During this year he published four papers in Annalen der Physik whose aftermath was prodigious.

The first scientific study of diffusion was performed by a Scottish chemist Thomas Graham (1805 -1869). His research on diffusion was conducted between 1828 and 1833. Here we quote his first paper: [ . . . ] the experimental information we possess on the subject amounts to little more than the well established fact, that gases of different nature, when brought into contact, do not arrange themselves according to their density, the heaviest undermost, and the lighter uppermost, but they spontaneously diffuse, mutually and equally, through each other, and so remain in the intimate state of mixture for any length of time.

In 1867 James Maxwell ‡ estimated the diffusion coefficient of CO 2 in the air using the results (measurements) of Graham. The resulting number was obtained within 5% accuracy. It is pretty impressive.

A. Fick § hold a chair of physiology in Würzburg for 31 years. His main contributions to Physics were made during a few years around 1855. When he was 26 years old he published a paper on diffusion establishing the now classical Fick's diffusion law. Fick did not realize that dissolution and diffusion processes result from the movement of separate entities of salt and water. He deduced the quantitative law proceeding in analogy with the work of J. Fourier [START_REF] Fourier | Théorie analytique de la chaleur[END_REF] who modeled the heat conduction: [ . . . ] It was quite natural to suppose that this law for diffusion of a salt in its solvent must be identical with that according to which the diffusion of heat in a conducting body takes place; upon this law Fourier founded his celebrated theory of heat, and it is the same that Ohm applied [ . . . ] to the conduction of electricity [ . . . ] according to this law, the transfer of salt and water occurring in a unit of time between two elements of space filled with two different solutions of the same salt, must be, ceteris partibus, directly proportional to the difference of concentrations, and inversely proportional to the distance of the elements from one another. * Albert Einstein (1879 -1955) was a German theoretical Physicist. This personality does not need to be introduced. † Annus Mirabilis comes from Latin and stands for "extraordinary year" in English or "Wunderjahr" in German.

‡ James Clerk Maxwell (1831 -1879) was a Scottish Physicist, famous for Maxwell equations.

§ Adolf Fick was the author of the first treatise of Die medizinische Physik (Medical Physics) (1856) where he discussed biophysical problems, such as the mixing of air in the lungs, the work of the heart, the heat economy of the human body, the mechanics of muscular contraction, the hydrodynamics of blood circulation, etc.

Going along this analogy, Fick assumed that the flux of matter is proportional to its concentration gradient with a proportionality factor κ, which he called "a constant dependent upon the nature of the substances". Actually, Fick made an error in (minus) sign which introduces the anti-diffusion and leads to an ill-posed problem. This error was not a source of difficulty for Fick since he analyzed only steady states in accordance with available experimental conditions at that time.

The fundamental article [START_REF] Einstein | Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen[END_REF] devoted to the study of Brownian * motion and entitled "On the Motion of Small Particles Suspended in a Stationary Liquid, as Required by the Molecular Kinetic Theory of Heat" by Einstein was published on 18 th July 1905 in Annalen der Physik. By the way, it is the most cited Einstein's paper among four works published during his Annus Mirabilis. This manuscript is of capital interest to us as well. Let us quote a paragraph from [START_REF] Einstein | Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen[END_REF]:

In this paper it will be shown that, according to the molecular kinetic theory of heat, bodies of a microscopically visible size suspended in liquids must, as a result of thermal molecular motions, perform motions of such magnitudes that they can be easily observed with a microscope. It is possible that the motions to be discussed here are identical with so-called Brownian molecular motion; however, the data available to me on the latter are so imprecise that I could not form a judgment on the question [ . . . ]

Einstein was the first to understand that the main quantity of interest is the mean square displacement X 2 (t) and not the average velocity of particles. Taking into account the discontinuous nature of particle trajectories, the velocity is meaningless.

Before publication of [START_REF] Einstein | Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen[END_REF] atoms in Physics were considered as a useful, but purely theoretical concept. Their reality was seriously debated. W. Ostwald † , one of the leaders of the anti-atom school, later told A. Sommerfeld ‡ that he had been convinced of the existence of atoms by Einstein's complete explanation of Brownian motion.

Without knowing it, Einstein in answered a question posed the same year by K. Pearson § in a Letter published in Nature [START_REF] Pearson | The Problem of the Random Walk[END_REF]: Can any of your readers refer me to a work wherein I should find a solution of the following problem, or failing the knowledge of any existing solution provide me with an original one? I should be extremely grateful for aid in the matter.

A man starts from a point O and walks ℓ yards in a straight line; he then turns through any angle whatever and walks another ℓ yards in a second straight line. He repeats this process n times. I require the probability that * Robert Brown (1773 -1858) was a Scottish Botanist. † Friedrich Wilhelm Ostwald (1853 -1932) was a Latvian Chemist. He received a Nobel prize in Chemistry in 1909 for his works on catalysis.

‡ Arnold Johannes Wilhelm Sommerfeld (1868 -1951) was a German theoretical Physicist. He served as the PhD advisor for several Nobel prize winners.

§ Karl Pearson (1857 -1936) was an English Statistician.

after these n stretches he is at a distance between r and r + dr from his starting point, O.

The problem is one of considerable interest, but I have only succeeded in obtaining an integrated solution for two stretches. I think, however, that a solution ought to be found, if only in the form of a series in powers of 1 n , when n is large.

The expression "random walk " was probably coined in Pearson's Letter [START_REF] Pearson | The Problem of the Random Walk[END_REF].

Let us follow Einstein's study of the Brownian motion. Consider a long thin tube filled with water (it allows us to consider only one spatial dimension). At the initial time t = 0 we inject a unit amount of ink at x = 0. Let u(x, t) denote the density of ink particles in location x ∈ R and at time t 0. So, initially we have

u(x, 0) = δ(x) , (C.1)
where δ(x) is the Dirac distribution centered at zero. Then, suppose that the probability density of the event that an ink particle moves from x to x + ∆x in a small time ∆t is ρ(∆x, ∆t). Then, we have

u (x, t + ∆t) = ˆ+∞ -∞ u(x -∆x, t) ρ(∆x, ∆t) d ∆x .
Assuming that solution u(x, t) is smooth, we can apply the Taylor * formula:

u (x, t + ∆t) = ˆ+∞ -∞ u -u x • ∆x + 1 2 u xx • ∆x 2 + • • • ρ(∆x, ∆t) d ∆x . (C.2)
Now let us recall that ρ is a Probability Density Function (PDF). Thus, Normalisation: ´+∞ -∞ ρ(∆x, ∆t) d ∆x = 1 Symmetry: we can assume that the diffusion process is symmetric in space, i.e. ρ(∆x, ∆t) = ρ(-∆x, ∆t) , ∀∆x, ∆t 0 .

Consequently,

ˆ+∞

-∞ ∆x ρ(∆x, ∆t) d ∆x = 0 .

Variance: We can also assume that the variance if finite and Einstein assumed additionally that the variance is linear in ∆t:

ˆ+∞ -∞ ∆x 2 ρ(∆x, ∆t) d ∆x = D∆t ,
where D > 0 is the so-called diffusion constant. * Brook Taylor (1685 -1731) was an English Mathematician.

By incorporating these results into equation (C.2) and rearranging the terms, it becomes

u(x, t + ∆t) -u(x, t) ∆t = 1 2 D u xx (x, t) + • • • .
Taking the limit ∆t → 0 we obtain straightforwardly

u t = 1 2 D u xx . (C.3)
The Initial Value Problem (C.1) for this linear parabolic equation (C.3) can be solved exactly:

u(x, t) = 1 √ 2πDt e -x 2 2Dt .
The last solution is also known as the Green function. However, the main result of Einstein's paper [START_REF] Einstein | Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen[END_REF] is the following formula:

D = RT f N a , (C.4)
where R: the ideal gas constant, i.e. R ≈ 8.3144598 J K • mol T : the absolute temperature f : the friction coefficient N a : Avogadro's * number, i.e. N a ≈ 6.022140857 mol -1 Formula (C.4) along with the observation of the Brownian motion enabled J. Perrin † to produce the first historical estimation of Avogadro's constant N a ≈ 6 mol -1 .

Exercise 1. Read a fascinating paper by A. Einstein on the Hydrodynamics of tea leaves [START_REF] Einstein | Die Ursache der Mäanderbildung der Flußläufe und des sogenannten Baerschen Gesetzes[END_REF]. By the way, the Author attests that Tea is very stimulating for intellectual activities.

D. Monte-Carlo approach to the diffusion simulation

From considering the history and physical modelling of Brownian motion we naturally come to the numerical simulation of parabolic PDEs using stochastic processes. It falls * Lorenzo Romano Amedeo Carlo Avogadro di Quaregna e di Cerreto (1776 -1856) was an Italian scientist.

† Jean Baptiste Perrin (1870 -1942) was a French Physicist who was honoured with the Nobel prize in Physics in 1926 for the confirmation of the atomic nature of matter (through the observation of Brownian motion).

where V (x) is a function representing the amount of external cooling (if V (x) 0, and external heating if V (x) < 0, but we do not consider this case) at point x. Then, we have the following Theorem 3 (Feynman-Kac formula). Let V (x) be a non-negative continuous function and let u 0 (x) be bounded and continuous. Suppose that u(x, t) is a bounded function that satisfies equation (D.1) along with the initial condition

u(x, 0) = u 0 (x) , (D.2) then, u(x, t) = E exp - ˆt 0 V W s ds u 0 W t , (D.3)
where {W t } t 0 is a Brownian motion starting at x.

We have to define rigorously also what the Brownian motion is:

Definition 1 (Brownian motion). A real-valued stochastic process t → W(t) (or a random curve) is called a Brownian motion (or Wiener process) if (1) Almost surely W(0) = 0 (2) W(t) -W(s) ∼ N (0, t -s) ≡ √ t -s N (0, 1), ∀t s 0 (3)
For any instances of time v > u > t > s 0 the increments W(v) -W(u) and W(t) -W(s) are independent random variables.

In particular, from point (2) by choosing s = 0 it follows directly that E W(t) = 0, E W 2 (t) = t . More informally we can say that the Brownian motion is a continuous random curve with the largest possible amount of randomness.

The Theorem above can be proved even under more general assumptions, but the formulation given above suffices for most practical situations. For instance, the functions V (x) and u 0 (x) may have isolated discontinuities and the Feynman-Kac formula will be still valid. Another interesting corollary of Feynman-Kac's formula is the uniqueness of solutions to the Initial Value Problem (IVP) (D.1): Corollary 1. Under the assumptions of Theorem 3, there is at most one solution to the heat equation (D.1), which satisfies the initial condition (D.2). Namely, it is given by Feynman-Kac formula (D.3).

The main (computational) advantage of Feynman-Kac formula is that it can be straightforwardly generalized to the arbitrary dimension:

Theorem 4 (Feynman-Kac formula in d dimensions). Let V : R d → [0, +∞) and u 0 : R d → R be continuous functions with bounded initial condition u 0 (x). Suppose that u(x, t) is also a bounded function, which satisfies the following Partial Differential Equation (PDE):

u t = 1 2 d i = 1 ∂ 2 u ∂x 2 i -V (x) • u ,
and the initial condition u(x, 0) = u 0 (x) . Then, u(x, t) = E exp -ˆt 0 V W s ds u 0 W t , where {W t } t 0 is a d-dimensional Brownian motion starting at x.

We give here another useful generalization of the Feynman-Kac method. For the sake of simplicity we return to the one-dimensional case. Consider a stochastic process X t , which satisfies the following Stochastic Differential Equation (SDE): dX t = α(X t ) dt + σ(X t ) dW t , X 0 = x , (D.4)

where α(x) is the local drift and σ(x) is called the local volatility in financial applications.

We assume here that α(x) and σ(x) are globally Lipschitz continuous and grow linearly in space at most. Then, the function given by formula (D.3) satisfies the generalized diffusion equation

u t = α(x) u x + 1 2 σ 2 (x) u xx -V (x) u(x) ,
together with the initial condition (D.2). So, the resulting numerical algorithm is very simple:

(1) Generate N trajectories of the Brownian motion

{W k s } 1 k N s ∈ [0, t]
(2) Compute N solutions to the SDE (D.4) using the Euler-Maruyama * method [START_REF] Higham | An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations[END_REF], for example (3) Compute the solution value using representation (D.3). The mathematical expectation E [•] is replaced by the simple arithmetic average according to the Monte-Carlo approach.

Remark 10. Notice that the step (2) can be omitted if we solve the linear parabolic equation without the drift term and with constant diffusion coefficient. The simplest version of the Feynman-Kac method for the simplest initial value problem

u t = 1 2 u xx , u(x, 0) = u 0 (x) , looks like u(x, t) ≈ 1 N N n = 1 u 0 (ξ n ) , ∀n : ξ n ∼ N (x, √ t) .
It is probably the simplest way to estimate solution value in a given point (x, t).

The Feynman-Kac method inherits all advantages (and disadvantages) of Monte-Carlo methods. Comparing to grid-based methods it enjoys also the locality property. Namely, if you want to compute your solution in a given point (x, t), you do not need to know the solution in neighbouring sites (locations). So, if you want to compute numerically the solution only in a few specific points, please, consider this method. It can be competitive with more conventional approaches. Moreover, the financial industry already appreciated the power of these methods.

D.1. Brownian motion generation

In this Section we provide a simple Matlab code to generate M sample Brownian motion realizations. It can be used as a building block for SDE solvers and practical implementations of the Feynman-Kac method described above: A sample output result (with precisely the same parameters) of this script is depicted in Figure 10.

It is quite straightforward to employ the same Matlab script to generate Brownian paths in d dimensions (just take M = m × d, where m ∈ N and then regroup matrix W in m d-dimensional trajectories). A few realizations of the Brownian motion in two spatial dimensions are depicted in Figure 11.

E. An exact non-periodic solution to the 1D heat equation

In this Appendix we provide an exact solution (see [START_REF] Fornberg | A practical guide to pseudospectral methods[END_REF]Section §7.2]) to the (dimensionless) heat equation Remark 11. Equation (E.1) (along with boundary conditions) models the temperature variation in soil under periodic boundary forcing (modeling seasonal temperature variations). Solution (E.2) explains also that for many types of soils, there is a phase shift of seasonal temperature at certain depths (warmest in winter and coolest in summer). See Figure 13 for an illustration. It provides a theoretical explanation why a cave works in practice.

F. Some popular numerical schemes for ODEs

The discussion of numerical methods for PDEs cannot be complete if we do not provide some basic techniques for the time marching. Basically, if we follow the Method Of Lines (MOL) [START_REF] Kreiss | Method of lines for hyperbolic equations[END_REF][START_REF] Reddy | Stability of the method of lines[END_REF][START_REF] Schiesser | Method of lines solution of the Korteweg-de vries equation[END_REF][START_REF] Shampine | ODE solvers and the method of lines[END_REF], the PDE is discretized first in space, then the resulting system of coupled ODEs has to be solved numerically. This document would not be complete if we did not provide any indications on how to do it. Consider for simplicity an Initial Value where the dot over a function denotes the derivative with respect to time, i.e. u ≡ du dt .

There is a number of time marching schemes proposed in the literature. We refer to [START_REF] Hairer | Geometric Numerical Integration[END_REF][START_REF] Hairer | Solving ordinary differential equations: Nonstiff problems[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems[END_REF] has actually infinitely many solutions. Two examples are u(t) ≡ 0 and u(t) = t 2 4 . Obviously, the right-hand side does not satisfy the Lipschitz condition which guarantees the uniqueness.

Figure 1 .

 1 Figure 1. Periodisation of a smooth continuous function defined on [-π, π].

Figure 2 .

 2 Figure 2. Interpolation of the Runge function f (x) = 1 1 + 25x 2 on a sequence of successively refined of uniform grids. One can observe the divergence phenomenon close to the interval boundaries x = ±1.

Figure 3 .

 3 Figure 3. Interpolation of the Runge function f (x) = 1 1 + 25x 2 on a sequence of successively refined of Tchebyshev's grids. One can see that the oscillations present in Figure 2 disappear and the interpolant seems to converge to the interpolated function f (x) .

Figure 4 .

 4 Figure 4. Illustration of the aliasing phenomenon: two Fourier modes are indistinguishable on the discrete grid. The modes represented here are cos(x) and cos(9x) and the discrete grid is composed of N = 11 equispaced points on the segment [-π, π].

1

  function w_hat = AntiAlias(u_hat, v_hat) N/2; % 3/2th rule 4 u_hat_pad = [u_hat(1:N/2) zeros(1, M-N) u_hat(N/2+1:end)]; 5 v_hat_pad = [v_hat(1:N/2) zeros(1, M-N) v_hat(N/2+1:end)]; 2 * [w_pad_hat(1:N/2) w_pad_hat(M-N/2+1:M)]; 11 end % AntiAlias()

Figure 5 .

 5 Figure 5. Comparison of various numerical approaches to determine the spectral expansion coefficients for a linear BVP (3.1).

1 l = 1 Figure 6 .

 16 Figure 6. Spectral solution to the linear heat equation (4.1) at t = 5.0 , with ν = 10 -2 and the initial condition is u 0 (x) = sech 2 (10 x) .

  the solution at final time: 16 uT = real(ifft(exp(-nu * k2 * T). * u0_hat));

1 l = 1 5 6 8 9

 158 .0; % half-length of the domain 2 N = 128; % number of Fourier modes 3 dx = 2 * l/N; % distance between two collocation points 4 x = (1-N/2:N/2) * dx; % physical space discretization dk = pi/l; 7 k = [0:N/2 1-N/2:-1] * dk; % vector of wavenumbers arg = pi * (x + 1); 10 sar = sin(arg); 11 car = cos(arg); car2 = car. * car; 12 esa = exp(sar); 13 u = sar. * esa; 14 uhat = fft(u); 15 16 % numerical derivatives: 17 up = real(ifft(1i * k. * uhat)); 18 upp = real(ifft(-k.^2. * uhat)); 19 uppp = real(ifft(-1i * k.^3. * uhat)); 20 21 % exact derivatives: 22 pi2 = pi * pi; pi3 = pi2 * pi; 23 u1 = pi * car. * (1 + sar). * esa; 24 u2 = pi2 * (car2. * (sar + 3) -1 -sar). * esa; 25 u3 = pi3 * car. * (car2. * (sar + 6) -4 -7 * sar). * esa;

Figure 7 .

 7 Figure 7. Numerical differentiation of a periodic function (4.9) using the

•

  The following book by N. Trefethen * has two major advantages: (i) conciseness and (ii) collection of Matlab programs which comes along. I would recommend it as the first reading on pseudo-spectral methods. At least you will learn how to program efficiently in Matlab (or in Octave, Scilab, etc. )[START_REF] Trefethen | Spectral methods in MatLab[END_REF]:-Trefethen, L. N.(2000). Spectral methods in MatLab. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA. * Lloyd Nick Trefethen (1955 -20..), an American/British numerical analyst.

2 P

 2

Figure 8 .

 8 Figure 8. The first five Tchebyshev polynomials. Notice that T n (-1) = (-1) n and T n (1) ≡ 1.

Figure 9 .

 9 Figure 9. Collocation Trefftz method using the Green function G(P, Q):{P i } 13i=1 ⊆ R d are collocation points and {Q i } 13 i=1 ⊆ R d are the sources put outside of the domain Ω in order to avoid singularities. Then, the approximate solution is sought as a linear combination of functions G(P, Q i ) 13 i=1 . The unknown expansion coefficients are found so that to satisfy the boundary conditions exactly in collocation points {P i } 13 i=1 . The outer points {Q i } 13 i=1 are not needed if we know a T-complete set of functions for the operator L .

1

  M = 100; % number of paths 2 N = 1000; % number of steps 3 T = 1; % final simulation time 4 dt = T/N; % time step 5 dW = sqrt(dt) * randn(M, N); 6 W = cumsum(dW, 2);

u t - 1 9π 2 uFigure 11 .

 211 Figure 11. A few Brownian paths in two spatial dimensions. The red circle depicts the starting point (0, 0).

Figure 12 .

 12 Figure 12. Space-time plot of the solution (E.2) to the linear heat equation (E.1). The time window is t ∈ [0, 9π].

Figure 13 .

 13 Figure 13. Solution (E.2) shown at the surface x = 0 and at depth x = 3 4 . The goal is to illustrate the phase shift between two curves.

  2]: -Boyd, J. P. (2000). Chebyshev and Fourier Spectral Methods. (Dover Publications, New York) (2 nd Ed.). • This book is probably my favourite one. It represents a good balance between the theory and practice and this book arose as an extended version of a previously published review paper. The present lecture notes were inspired in part on this book as well [11]: -Fornberg † , B. (1996). A practical guide to pseudospectral methods. Cambridge: Cambridge University Press. • The Author discovered this book only recently. So, he has still to read it, but from the first sight I can already recommend it:

	In general, the Author's suggestions (according to his personal taste and vision) for
	scientific literature and software are collected in a single document, which is continuously
	expanded and completed:
	https://github.com/dutykh/libs/
	Good luck and have a nice reading!
	⋆ ⋆ ⋆

-Peyret, R.

(2002)

. Spectral Methods for Incompressible Viscous Flow. Springer-Verlag New York Inc.

• Finally, I discovered also these very clear and instructive Lecture notes

[START_REF] Uecker | A short ad hoc introduction to spectral methods for parabolic PDE and the Navier-Stokes equations[END_REF]

. Lectures were delivered in 2009 at the International Summer School 'Modern Computational Science' in Oldenburg, Germany. They can be freely downloaded at the following URL address:

http://www.staff.uni-oldenburg.de/hannes.uecker/pre/030-mcs-hu.pdf * John Boyd (19.. -20..), an American numerical analyst, meteorologist and occasional science fiction writer.

† Bengt Fornberg (19.. -20..), a Swedish/American numerical analyst.

* We can apply the Pythagoras theorem since the aliasing error R N [u] contains the Fourier modes with numbers | k | m, while the reminder u -T N [u] contains only the modes with | k | > m. Thus, they are orthogonal.

* The advantage of the collocation approach for nonlinear (and variable coefficients) problems becomes even more flagrant.† With a little advantage towards the collocation method. However, this conclusion is not general at all. It is based only on this particular problem.

Acknowledgments

The Author would like to thank Professor Didier Clamond (University of Nice Sophia Antipolis, France) for introducing me to Fourier-type pseudo-spectral methods. I would like also to thank Professor Laurent Gosse (CNR, Italy) who brought Author's attention to Trefftz methods. Special thanks go to my colleague Prof. Paul-Éric Chaudru De Raynal (LAMA, Université Savoie Mont Blanc, France) for his help with Monte-Carlo simulations. Finally, the Author would like to thank Professors Nathan Mendes and Marx Chhay for giving me an opportunity to deliver these lectures in Brazil.

into the large class of Monte-Carlo and quasi-Monte-Carlo methods [START_REF] Caflisch | Monte Carlo and quasi-Monte Carlo methods[END_REF]. They are based on the so-called Law of large numbers, which can be informally stated * as

where ξ n are independent, identically distributed random variables and g(•) is a real-valued continuous function. Here E [•] stands for the mathematical expectation and the convergence is understood in the sense convergence in probability distribution † . The proof is based on the Tchebyshev inequality in Probabilities, but we do not enter into these details. The main advantages of Monte-Carlo methods are

• Simplicity of implementation • Independent of the problem dimension (they do not suffer of the curse of dimensionality) On the other hand, the convergence rate is given by the Central Limit Theorem [START_REF] Jaynes | Probability Theory[END_REF] and it is rather slow, i.e. O(N -1/2 ), even if some acceleration is possible thanks to some adaptive procedures [START_REF] Lapeyre | A framework for adaptive Monte-Carlo procedures[END_REF] (such as low-discrepancy sequences, variance reduction and multi-level methods). The large deviation theory guarantees that the probability of falling out of a fixed tolerance interval decays exponentially fast.

To the Author knowledge, the method we are going to describe below was first proposed by R. Feynman ‡ in order to solve numerically the linear Schrödinger § equation in 1940s. In fact, he noticed that the Schrödinger equation can be solved by a kind of average over trajectories. This observation led him to a far-reaching reformulation of the quantum theory in terms of path integrals. Upon learning Feynman's ideas, M. Kac ¶ (Feynman's colleague at Cornell University) understood that a similar method can work for the heat equation as well. Later it became Feynman-Kac method. Unfortunately, this method is not implemented in any commercial software (for PDEs, the financial industry is using these methods since many decades) and it is not described in classical textbooks on PDEs and in the books on numerical methods for PDEs.

Consider first the following simple heat equation:

) * We choose this particular form, since it is suitable for our exposition. † A sequence of random variables {ξ n } ∞ n=1 converges in probability distribution towards the random variable ξ if for ∀ ε > 0 lim

This fact can be denoted as ξ n P ξ. ‡ Richard Phillips Feynman (1918 -1988) (E.1) described above for ∀t > 0 is given by *

This analytical solution can be used as the reference solution in order to validate your non-periodic numerical codes. Notice that the second term u Σ (x, t) (i.e. the infinite series) vanishes uniformly in space, i.e.

This term u Σ is needed to enforce the initial condition. Consequently, the long time behaviour of the solution (E.2) is given by u (x, t).

as exhaustive references on this topic. Below we provide some most popular (subjectively) schemes.

Forward Euler: (explicit, first order accurate)

Backward Euler: (implicit, first order accurate)

Adams * -Bashforth-2: (explicit, second order accurate)

Adams-Bashforth † -3: (explicit, third order accurate)

Adams-Moulton ‡ -1 or the trapezoidal rule: (implicit, second order accurate)

Adams-Moulton-2: (implicit, third order accurate)

. * John Couch Adams (1819 -1892), a British Mathematician and Astronomer. He predicted the existence and position of the planet Neptune.

† Francis Bashforth (1819 -1912) is a British Applied Mathematician working in the field of Ballistics. However, his famous numerical scheme was proposed in collaboration with J. C. Adams to study the drop formation.

‡ Forest Ray Moulton (1872 -1952) is an American Astronomer. There is a crater on the Moon named after him.

* Martin Wilhelm Kutta (1867 -1944) is a German Mathematician who worked in the field of Fluid Mechanics and Numerical Analysis.

Two-stage Runge-Kutta * : (explicit, second order accurate)

• α = 1 2 : mid-point scheme

Runge-Kutta-4 (RK4): (explicit, fourth order accurate)

We do not discuss here the questions related to the stability of these schemes. This topic is of uttermost importance but out of scope of these Lecture notes. Moreover, we skip also the adaptive embedded Runge-Kutta schemes which can choose the time step to meet some prescribed accuracy requirements [7].

In order to illustrate the usage of RK4 scheme (sometimes called the Runge-Kutta scheme) we take a simple nonlinear logistic equation:

It is not difficult to check that the exact solution to this IVP is

The Matlab code used to study numerically the convergence of the RK4 scheme on the nonlinear equation (F.1) is given below: † Karl Heun (1859 -1929) is a German Mathematician. The numerical result is shown in Figure 14. One can clearly observe the 4 th order convergence of the RK4 scheme applied to a nonlinear example. The convergence is broken only by the rounding effects inherent to the floating point arithmetics.

F.1. Existence and unicity of solutions

Normally, before attempting to solve a differential equation, one has to be sure that the mathematical problem is well-posed. The mathematical notion of the well-posedness was proposed by J. Hadamard * [START_REF] Hadamard | Sur les problèmes aux dérivées partielles et leur signification physique[END_REF] and it includes three points to be checked:

(1) Existence ( 2

) Uniqueness (3) Continuous dependence on the initial condition and other problem parameters

The last point is usually more difficult to be proven theoretically. However, for 3D Navier † -Stokes ‡ equations even the global existence and smoothness of solutions (i.e.the first point) poses already a Millennium problem: http://www.claymath.org/millennium-problems/navier-stokes-equation N, points Below we give some theoretical results which address the well-posedness conditions for an Initial Value Problem (IVP) for a scalar equation:

The generalization to systems of Ordinary Differential Equations (ODEs) is straightforward.

Theorem 5 (Existence). If function f (t, u) is continuous and bounded in a domain (t, u) ∈ D ⊆ R t × R u , then through every point (t 0 , u 0 ) ∈ D passes at least one integral curve of equation (F.2).

Theorem 6 (Prolongation). Let function f (t, u) is defined and continuous in domain D . If a solution u = φ(t) to problem (F.2) exists in the interval t ∈ [t 0 , α) cannot be prolongated beyond the point t = α, then it can happen only for one of three following reasons:

(

When t → α -0, distance between the point t, φ(t) to the boundary ∂D goes to zero.

Theorem 7 (Uniqueness, [START_REF] Osgood | Beweis der Existenz einer Lösung der Differentialgleichung dy dx = f (x, y) ohne Hinzunahme der Cauchy-Lipschitz'schen Bedingung[END_REF]). If function f (t, u) is defined in domain D and for any pair of points (t 1 , u 1 ) , (t 2 , u 2 ) ∈ D satisfies the condition

where ω(u) > 0 is continuous for 0 < u U and

then through any point (t 0 , u 0 ) ∈ D passes at most one integral curve of equation (F.2).

Suitable functions ω(u) which satisfy the conditions of the Theorem are

Above K is a positive constant. If we take ω(u) = K u then condition (F.3) becomes the well-known Lipschitz * condition for function f (t, u) in the second variable u . In order to satisfy the Lipschitz condition, it is sufficient for function f (t, u) to be defined in a domain D convex in u and to have a bounded derivative ∂f ∂u in this domain. Later, Wintner showed that Osgood † 's Theorem conditions are sufficient for the convergence of Picard ‡ 's iterations to a local solution on a sufficiently small interval [START_REF] Wintner | On the Convergence of Successive Approximations[END_REF].

F.1.1 Counterexamples

The Theorem above gives the existence of solutions only locally in time t 0 , t 0 + δ . The only reason for a solution not to exist beyond the given interval is the blow-up (i.e. the solution becomes unbounded). Otherwise, the solution exists globally. For instance, the following problem u = 1 + u 2 , u(0) = 0 , has the unique exact solution u(t) = tan(t) which exists only in the interval 0, † William Fogg Osgood (1864 -1943) was an American Mathematician born un Boston, MA. He studied in Universities of Göttingen and Erlangen, Germany. He made contributions to Mathematical and Complex Analysis, Conformal Mappings.

‡ Charles Émile Picard (1856 -1941) was a French Mathematician who made contributions to the Mathematical Analysis. He was married to Marie, a daughter of his Professor Charles Hermite.