Spreading of particles in some displacement probability field

SPREADING OF PARTICLES IN SOME

DISPLACEMENT PROBABILITY FIELD

T ^ ln\ K05

SPREADING OF PARTICLES IN SOME

DISPLACEMENT PROBABILITY FIELD

SPREADING OF PARTICLES IN SOME

DISPLACEMENT PROBABILITY FIELD

Proefschrift

ter verkrijging van de graad van doctor in de technische wetenschappen aan de Technische Universiteit Delft, op gezag van de

Rector Magnificus, prof.dr. J.M. Dirken,

in het openbaar te verdedigen ten overstaan van

het College van Dekanen op dinsdag 11 november 1986 te 14.00 uur

door

RAPHAEL SOERJADI

geboren te Muntilan (Indonesië) wiskundig ingenieur

Dit proefschrift is goedgekeurd door de promotor prof.dr.ir. A . J . Hermans

Dedicated to Floor and Bima Kresna -Wibi-Ardjoena-Wisnoe

S T E L L I N G E N

1. Als een wolk van deeltjes, op tijdstip t=0 gaussisch verdeeld met gemiddelde y en covariantiematrix C als parameters, wordt losgelaten in een n-dimensio-naal gaussisch verplaatsingskansveld met veldcon-stanten a en B, dan blijft de verdeling van de deel­ tjes gaussisch met op tijdstip t de gemiddelde u + at en covariantiematrix C + 2Bt.

■=o — o Dit proefschrift, paragraaf 11.2.

2. Als het verplaatsingskansveld in een n-dimensionale ruimte bestaat uit een rooster van M onafhankelijke soorten kanalen zodanig dat M £ n, dan is de daar­

bijbehorende tweede orde differentiaalvergelijking ontaard.

Dit proefschrift, paragraaf 14.1.

3. Als alle atomen van een differentie-molecule, ont­ wikkeld voor het numeriek oplossen van een partiële differentiaalvergelijking van het parabolisch-achtige type, groter dan of gelijk zijn aan nul, dan is het daarbijbehorende numerieke proces stabiel. Behalve voor de tweede orde vergelijking, is deze stelling in het algemeen niet omkeerbaar.

Dit proefschrift, paragraaf 16.4.

4. De momenten M van orde p+q (met p,q = 0,1,2,...) van een willekeurige veelhoek (met hoekpunten 1,2,..)

>=

ff xP y^ dx dy

kunnen worden berekend als de algebraische som van de momenten van de driehoeken gevormd door de oorsprong en een zijde van de veelhoek. Voor de driehoek met zijde van hoekpunt 1 naar hoekpunt 2 wordt gevonden :

p! q! M (pq)_ 12 (p + q + 2) ! xl X2 i+j\ /P+q-i-J1 p - i 1 p_ 1 J xl x2 yl y2

Dit biedt de mogelijkheid om de oppervlakte-inte­ graal met behulp van een computer uit te rekenen door slechts de contour (dus de hoekpunten) van de veelhoek te geven.

*R.Soerjadi, HERON, 1966, No.3; J.II.J.Almering, IIERON, 1967, No.1; R.Soerjadi, HERON, 1968, No.5.

5. Mechanica kan worden gedefinieerd als de leer der verplaatsingen van deeltjes-verzamelingen. Toege­ paste Mechanica (Engels: Structural Mechanics) is dat deel van de mechanica waarbij de verplaatsingen (al of niet verwaarloosbaar) klein zijn.

6. Goed doceren impliceert goed doseren.

7. Het oude Javaanse spreekwoord: O /

(Ut A&MiKn ai i&inn&tui*

béja, waarbij de e als in 'we' en de a als de o in 'bok' wordt gelezen, en vrij vertaald: Voor wel­ vaart/welzijn moet je betalen), dat vermoedelijk tot stand is gekomen in de economische en geeste­ lijke betekenis, krijgt, naar de huidige inzichten, ook een ecologisch karakter.

8. Het verdient aanbeveling om het wachttijdprobleem van patiënten aan te pakken. Dit zal velerlei voor­ delen opleveren, waaronder ook economische.

9. Het begrip 'absoluut gehoor' in de muziek is niet echt absoluut in de absolute betekenis van het woord. Immers, de toon a uit het ééngestreept oc­ taaf, bij voorbeeld, is tegenwoordig (Londen 1939) vastgesteld op 440 Hz, terwijl dat vroeger (bij­ voorbeeld: Parijs 1648) 374 Hz was. Dit betekent dat iemand met toen een absoluut gehoor nu de plank mis zou slaan.

10. In het schaakspel heeft de witspeler een licht voordeel ten opzichte van zijn tegenstander omdat hij mag beginnen. Een soortgelijk voordeel wordt in het (met schaken vergelijkbare) go-spel gecom­ penseerd door 5 'komi', dat zijn 5 gevangen stenen die in de eindstand worden verrekend.

Stelling: Het verdient aanbeveling om ook in het schaakspel, eerlijkheidshalve, een dergelijke com­ pensatie reglementair in te voeren, bij voorbeeld in de vorm van 1 a 2 % verschil in bedenktijd.

11. Sommige talen hebben lidwoorden (bijvoorbeeld het Nederlands) en andere niet (bijvoorbeeld de Bahasa Indonesia).

Stelling: Zuiver als informatie-verstrekkend com­ municatiemiddel zou de Nederlandse taal het heel goed zonder lidwoorden kunnen stellen. Dit is slechts een constatering en geen aanbeveling.

PREFACE

It was about a decade ago that my teacher at the Department of Civil Engineering, G. de Josselin de Jong, drew my attention to the problem of dispersion of particles. He himself has studied this topic,

resulting in many publications (among which: ref. 7 to 11). Among other things, he has developed a differen­ tial equation that describes the dispersion process (ref. 10). There was, however, some doubt about its correctness, for a reason that is not relevant here. At his request and in collaboration with him, the present study has been initiated. The aim of this study is mainly to develop the governing equation of the dispersion process, based upon the same model used by de Josselin de Jong, applying, however, a different mathematical method. So the present study has been based on the work of de Josselin de Jong, but arranged in another (for this purpose more con­ venient) way, and here and there filled up with some little ideas. The differential equation found in this study turns out to be exactly the same as the one developed by de Josselin de Jong. This strengthens my opinion that his conceptual idea was correct. I wish to express my gratitude to him for the papers he has allowed me to use, for the many personal communica­ tions and especially for the courage he gave me in tackling the problem.

A.J.Hermans from the Sub-department of Mathematics and Informatics has also helped and encouraged me in this study. I have got some inspirational ideas during

the many personal communications with him. I thank him for this.

I am also indebted to the Management of the Group of Structural Mechanics of the Department of Civil Engin­ eering (where I am working), in particular its chair­ man A.L.Bouma, for giving me the opportunity to carry out this study.

Besides, I am grateful to A.L.Bouma, H.W.Loof, and A. Verruijt, for their stimulating ideas and encourage­ ment.

C O N T E N T S

page

PREFACE 5

1. INTRODUCTION 11 2. PROBLEM FORMULATION AND GENERAL APPROACH 13

2.1. The main problem 13 2.2. General approach 13 3. PARTICLE DISTRIBUTION FUNCTION 15

4. MACRO-MODEL OF THE DISPLACEMENT PROBABILITY

FIELD 17 4.1. Setting up the macro-model 17

4.2. Properties of the displacement

probability field 22 5. ELEMENTARY VOLUME 25

5.1. Setting up the elementary volume 25 5.2. How to use the elementary volume 27 6. DEVELOPMENT OF THE GOVERNING INTEGRAL

EQUATION 30 7. EXAMPLE SATISFYING THE INTEGRAL EQUATION 36

7.1. One-dimensional 36 7.2. Two-dimensional 43 7.3. More-dimensional 53 8. DERIVATION TO DIFFERENTIAL EQUATION 63

9. REDUCTION TO SECOND-ORDER EQUATION 71 10. EXAMPLE SATISFYING THE SECOND-ORDER EQUATION 76

11. OTHER EXAMPLES 79 11.1. Instantaneous injection in a place

11.2. A Gaussian particle distribution as

the initial condition 84 11.3. Other simple initial conditions 89

11.4. A point injection during a time

interval 90 12. MICRO-MODEL OF THE DISPLACEMENT PROBABILITY

FIELD 92 12.1. Setting up the micro-model 92

12.2. Properties of the displacement

probability field 98 13. DEGENERATION OF SECOND-ORDER EQUATION 110

13.1. General approach 110

13.2. Example 1 112 13.3. Example 2 117 13.4. Example 3 119 14. (NON-)DEGENERATION DEDUCED FROM MICRO-MODEL 122

14.1. A theorem for some class of problems 122 14.2. Example leading to a degenerate

equation 127 14.3. Example governed by a non-degenerate

equation 130 15. NUMERICAL APPROACH (GENERAL REMARKS) 135

16. ONE-DIMENSIONAL PROBLEM TREATED BY

FINITE-DIFFERENCE METHOD 137 16.1. Second-order equation 137 16.2. Third-order equation 140 16.3. Fourth-order equation 152 16.4. Stability analysis 159 16.5. Worked example 170

16.6. Direct calculations of the moments 184

16.7. Review of the results 185

17. NUMERICAL SIMULATION 190 17.1. Introduction 190 17.2. One-dimensional 191 17.3. Two-dimensional 197 18. SUMMARIZING CONCLUSIONS 211 REFERENCES 214

1. INTRODUCTION

Spreading of particles appears in various problems in physics and engineering. Many physicists and engineers have studied this topic from specific problems with their specific physical properties, such as particles in ground water flow or in porous media. Lack of know­ ledge of these specific problems has forced me to study this phenomenon in a general way. The present study is an attempt to describe the spreading of particles mathematically, with a minimum of physics vaguely on the background.

On the one hand, this general approach has the advantage that the method can, with some adaption, perhaps be applied in many similar problems, such as diffusion, i.e. spreading of particles on molecular level, or perhaps even spreading of stars on galactic level. The latter is of course exaggerated, and is only in a manner of speaking. On the other hand, it has the dis­ advantage that it ignores the specific character of the physical problem, while this often forms the kernel of the problem. And as such it cannot be disregarded as it is an essential part of the problem.

The word 'dispersion' often has a specific meaning in some problems. In this study it has the same meaning as

'spreading of particles'. The latter is preferred, as it indicates the general character of the treatment. The fundamental concept of 'particles', generally used in studies of this kind, will also be applied here. Part­ icles are supposed to be small, whatever it means. Among other things it implies that an arbitrary number of

particles can be located at one point in space without disturbing each other. Something like saturation is out of the question. Besides, there are no mutual interac­ tions between the particles.

This study is a compound of preceding publications (ref. 14 to 2 3 ) . It has been based on the work of de Josselin de Jong (ref. 7 to 1 1 ) , but somewhat modified for the sake of convenience.

The first sections deal with the development of the governing equation. It turns out to be an integral equation, which can be reduced to a differential equa­ tion, if necessary. In the next sections, a simple exam­ ple is presented, for which an analytical solution of the integral or differential equation can be obtained. In case the problem cannot be solved analytically, a numerical approach can be considered. This study spends a section to the treatment of the finite difference method, adapted to the problem under consideration. An important question is how these (analytical and numerical) results can be verified. Strictly speaking, verification can only take place by 'reality', i.e. in a real process. Generally, however, we will do with an experiment. But, a physical experiment, even on scale in a laboratory, is expensive and difficult to carry out. Instead, a numerical simulation on a computer will

be executed. This numerical simulation can be considered as 'reality', the results of which can serve as a touch­

2. PROBLEM FORMULATION AND GENERAL APPROACH

2.1. The main problem

A cloud of particles, with a given initial distribution, has been released in space at a given time. This part­ icle cloud will be spread due to some scattering mech­ anism in that space. How is the distribution of the particle cloud at another time?

Apparently, this is an initial and boundary value prob­ lem. Instead of asking for the particle distribution at a particular time, we can question how particles pass through at a given place. This is another problem. In the first problem the particle distribution in space is asked at a given time, whereas the second problem is questioning the particle distribution in time at a given place. This second problem, however, although mathemat-, ically equivalent with the first one, is not very im­ portant practically. It will not be discussed in this study. For the sake of completeness, it will be glanced at occasionally.

2.2. General approach

Strictly speaking, the spreading of the particle cloud depends on the properties of the particles, on the structure of the space in which the particles are

released, and on eventual interactions between them. In this study the problem is simplified ultimately. All those difficult things are swept together and replaced by a model, a mathematical-physical model. This model

desVibes how particles are displaced. And, since these displacements have a random character, generally, the model can be called a displacement probability field. This means, roughly speaking, that a particle released in such a field will undergo a random displacement according to some probability function.

After having defined the model, the next thing to worry about is, which equation is governing the particle spreading process? Such an equation can be obtained by analysing the change of the particle distribution with respect to a suitable elementary volume. But, what is a suitable elementary volume? The answer depends on how the problem is posed. The main problem has the purpose to determine the particle distribution at a particular time. In this case we have to know how the particle distribution changes with time. Then we have to observe the particle distributions at two different times sep­ arated by a given time interval. So this time interval is an essential part of a suitable elementary volume. This will be treated in section 5. Another problem, for instance the second problem mentioned above, would need a different elementary volume.

From these observations, the governing equation can be derived. Then, the obtained equation should be solved, taking into account the given initial and boundary conditions. But before doing so, we first must define properly what is meant by particle distribution.

3. PARTICLE DISTRIBUTION FUNCTION

A particle cloud means that particles are distributed in space. This distribution can be denoted by a function of the place coordinates. But this is not enough, be­ cause the particle distribution is changing with time. So it is also a function of the time variable t. An­ other aspect is that the distribution function can be discreet or continuous. For instance, in one-dimensional space with coordinate x, the particle distribution can be denoted by the function U(x^»t) (i = 1,2,...) in the discreet case, or by u(x,t) in the continuous case. The total amount of particles is assumed to be constant. There are no sources or sinks in the field, and par­ ticles cannot be gained or lost anyhow. Then, without loss of generality, the total amount of particles can be normalized, i.e. made equal to unity. In other words, the particle distribution function can be defined such that

J2 UCxi.t) = 1 at any time t (3.1)

i

in the discreet case, or

oo

I u(x,t) dx = 1 at any time t (3.2)

— 0 0

in the continuous case. For the sake of simplicity, the integration domain of (3.2) is supposed to be the entire space. It may be replaced, if desired, by the given

re-gion. By this normalization, the functions U(x^.t) and u(x,t) can be interpreted as a probability function and a probability density function, respectively. In fact, they are conditional probability (density) functions.

For, they denote the probability for a particle of being located at x^ or x at a given time t. In probability

theory they should be written as UCx-j/t) and u(x/t), separating the two variables by a slash. But, in order to make clear that they are functions of both variables, place and time, the notations UCx.pt) and u(x,t) would have been preferred. On the other hand, however, the existence of a sum or integral of these functions at a particular place over a time domain is not necessary. Moreover, it has no physical meaning at all. Therefore we steer a middle course. These functions will be de­ noted by UCx.pt) and u(x;t), separating the two varia­ bles by a semicolon. This semicolon means that mathe­ matically the function is really, and thus can be treat­ ed as, a function of two independent variables, while physically it is a function of the first variable only, the second variable always having a particular value. The function u(x;t), generally, can be partly continu­ ous in the bounded space considered, i.e. continuous on some subsets of that space. This would require a more careful approach. For the sake of simplicity, we only consider either a continuous or a discreet case. Sometimes, even the discreet case will not be mentioned.

4. MACRO-MODEL OF THE DISPLACEMENT PROBABILITY FIELD

4.1. Setting up the macro-model

A particle, being released in space, will undergo a displacement. This displacement is a fundamental va­ riable, and can be either stochastic or deterministic. Another fundamental variable can be the time interval in which that displacement takes place. This time in­ terval too is either stochastic or deterministic. These two variables will be used to set up the model. Instead we could have chosen, for instance, the particle velo­ city as fundamental variable. This would lead to other models.

For the moment, we do not care how the displacement and the time interval have been built-up. In such a case, the model will be called a macro-model. If we cannot deduce directly the properties of the macro-model, then we are forced to investigate how the displacement has been built up, arriving at a so-called micro-model. This will be treated in section 12.

Once we have chosen the two fundamental variables above, we can distinguish the following four possibilities: (1) The particle undergoes a stochastic displacement

during a deterministic time interval.

(2) The particle undergoes a deterministic displacement during a stochastic time interval.

(3) The particle undergoes a deterministic displacement during an also deterministic time interval.

(4) The particle undergoes a stochastic displacement during an also stochastic time interval.

The first three possibilities can be understood easily. The fourth one, however, is somewhat more difficult. In this case there should be, necessarily, another deter­ ministic variable which can serve as parameter. Hereto we can use the dimensionless step magnitude. It is a step as well in space as in time. This step parameter connects the particle displacement and the time inter­ val in a similar way (see ref. 22). It will be discuss­ ed later on.

The four possibilities represent that many models of the field. It will turn out that the first model com­ plies directly with the main problem (see section 2 ) , in which the interest is focussed on the particle dis­ tribution in space at a particular time. The second model complies with the second problem, mentioned in section 2, in which the interest is focussed on the particle distribution in time at a particular place. The third model is entirely deterministic. In that case we are not dealing with a probability field anymore. Strictly speaking, it can be considered as a special case of the first or of the second model. It will not be discussed at all. The fourth model, on the contrary, can be regarded as a generalisation of the first or of the second one. It gives the ability of treating the place coordinates and the time variable in exactly the same manner, emphasizing their symme­ tric behaviour. Reversely, the first and second model can be considered as special cases of the fourth one. The fourth model can be used for solving the main prob­ lem as well as the second problem. In both cases the step parameter should be eliminated, such that their

final solutions do not contain the step parameter anymore.

For the sake of simplicity, consider the one-dimen­ sional space with Cartesian coordinate x. The particle displacement in x-direction can be denoted by X. The time variable has been denoted by t, then similarly the time interval in which X takes place, can be in­ dicated by T. In the first model, the variable X is stochastic, while T is deterministic. This can be made clear by figure 1. The other models are illustrated by the figures 2, 3, 4a, 4b. The arrows in the figures suggest that the particle paths are straight lines. This is not necessarily true, the paths can be tortu­ ous. But it does not count in this macro-model. Only the starting and arriving points are important.

t

X

x '

; X

0

*

Fig.I. A particle, released at point xQ at time tQ,

will undergo a stochastic displacement X during a deterministic time interval T.

t +T o

Fig. 2. A particle, released at point x at time t , will undergo a deterministic displacement X during a stochastic time interval T.

Fig. 3. A particle, released at point x at time t , will undergo a deterministic displacement X during a deterministic time interval T,

I »■

V-

til

Fig. 4a. A particle, released at point x at time t , will undergo a stochastic displacement X during a stochastic time interval T.

Po+P

Fig. 4b. A particle, released at point x after having made a step magnitude p , will undergo a stochastic displacement X in a deterministic step magnitude P.

4.2.

Properties of the displacement probability field

Now consider the first model in one-dimensional space. The properties of the displacement probability field can be represented by a probability function. This is a function of the stochastic displacement X and the de­ terministic time interval T in which the displacement takes place. Generally, the displacement X is depend­ ing on where and when the particle is released. So X is a function of x and t. For the sake of simplicity, how­ ever, this study is dealing with probability fields which are place and time independent. This means physically, that the spreading mechanism is covering homogeneously the x-space and does not change with the time t. If required, the method presented here could be extended to more complicated cases.

We are always dealing with conditional probability func­ tions, as T is deterministic. So F(X;T) means the con­ ditional probability for a particle, whereever and when­ ever it is released, of undergoing a displacement X during a given time interval T. The use of a semicolon has been explained before in section 3. It denotes ma­ thematically that F(X;T) is, and thus can be treated as, a function of X and T, but physically, it is a function of X only, with T having some particular value. The same reasoning applies for the conditional density

function f(X;T), where the displacement is continuous. The second model can be treated similarly. The corres­ ponding conditional probability density function then will be f(T;X). Consequently, the (normalized) particle distribution function should be denoted by u(t;x).

Ma-thematically, we only have to interchange the variables X and T, and the variables x and t. But physically, it has, of course, a different interpretation.

The third model is entirely deterministic. Formally, we can say that f(X;T) equals one for some value X0,

and equals zero for other values of X. Or, we may also say that f(T;X) equals one for some particular value TQ, and equals zero for other values of T. This model,

however, will not be discussed in this study.

The fourth model deals with the conditional probability density function f(X;P), according to figure 4b, or with f(X,T;P), according to figure 4a. The determinis­ tic step parameter P now plays the role of the time interval T, while T itself can be considered as a 'dis­ placement' in time. So this model can be treated as a more-dimensional first model, provided that the step parameter will be eliminated afterwards.

The story above presumes a given conditional probab­ ility (density) function. Generally, however, this conditional probability (density) function is not given, at least not given explicitly. A given expres­ sion of the displacement will do, although it will not be that easy. In such a case we can deduce its moments, if they exist, which determine implicitly the con­

ditional probability (density) function. More difficult is the case in which neither the expression of the dis­ placement nor the conditional probability (density) function is given. In that case, a close-up of the dis­ placement probability field will be necessary. It has been mentioned before that then we are forced to con­ sider the so-called micro-system. This micro-system

needs its own mathematical-physical model (see section 12). From this micro-model the properties of the dis­ placement probability field can be deduced.

5. ELEMENTARY VOLUME

5.1. Setting up the elementary volume

For the sake of simplicity, consider the first model in the one-dimensional space. Other cases can be treat­ ed similarly. Let the conditional probability for a particle of undergoing a displacement X-^ (i = 1,2,...) in a given time interval T be denoted by FCXj^T). The domain D(X;T) is the set of all possible values of X^ in the given time interval T.

D(X;T) = {X^x = {Xlt X2 *T (5.1.1)

Then, since a particle always undergoes some displace­ ment XJL in that given time interval T,

E F ( X±; T ) = 1 (5.1.2)

i

If X is continuous, i.e. X is an element of a continu­ ous set D(X;T), and the corresponding conditional prob­ ability density function is denoted by f(X;T), then (5.1.2) becomes

ƒf(X;T) dX = 1 (5.1.3) D(X;T)

Let Xm„„ and X_. be the maximum and minimum value of

IaclX mill

the set D(X;T), respectively. Then, in fact, we are considering an elementary volume of size (Xmax- Xm^n) .

This size depends on the given time interval T. It can also be said, that we are considering an elementary volume of size (Xraax- Xm i n, T) in the xt-space (see

figure 5 ) .

t

V

1

t

X -X '

max min

*min

....L-2&

X

max ^

X

T

As mentioned before in section 2, this elementary volume is suitable for treating the main problem. It gives the ability of observing the particle distribu­ tion in space at two different times separated by a given time interval T. Principally, an elementary volume should be small. And generally, it is. But exceptionally, in some concepts, it can be large, even infinite. Nevertheless, it still has to satisfy the weaker condition

lim X m 0

T-K)

(5.1.4)

Similar formulae apply for the other models and in more-dimensional cases.

5.2. How to use the elementary volume

It seems to be peculiar talking about this, but it is not. Because, strictly speaking, there are, in the model mentioned above, two different ways of using the elementary volume. As a matter of fact, figure 5 is suggesting the use of the elementary volume in a par­ ticular way, the first way, say. This is, by consider­ ing a particle at place xQ and at time tQ, and asking

where that particle is going to after a given time in­ terval T. The other way, to be called the second way, is by questioning where a particle at time t0 at a

given place x^ could have come from, a given time in­ terval T before. These are the two ways, in which the elementary volume can be used. It is illustrated by figure 6.

t

t +T

o

' o

max min

min

1 \ A

—L-AC-X

max

X

Fig. 6. Two ways of using the elementary volume: the 'going to' (6.1) and the 'coming from'(6.2) way (see next page).

V

T

max

min

JZEI

x. x

t

v

T

X

""I Zfö~~* 1

\ss

!

V

Fig. 6. Two ways of using the elementary volume: the 'going to' (6.1) (see previous page), and the

In both cases, the elementary volumes are of the same size, although they have been built up differently. The essential difference is denoted by the arrows in the figure. It can be characterized as the 'going to' and the 'coming from' way. It should be noted that in the third model these two ways are identical, because it is entirely deterministic. We might overlook these two possibilities, as we are concerned with deterministic problems, usually.

Principally, it does not matter, whether the state at time tQ will be compared with the one at time tQ+T, or

the state at time tQ-T with the one at time tQ. It will

turn out that only the second way can be used success­ fully in our problem. In the second way, we may take either tQ and tQ+T (see figure 6.2b) or tQ-T and tQ

6. DEVELOPMENT OF THE GOVERNING INTEGRAL EQUATION

Consider, for the moment, the one-dimensional case. Extension to more dimensions will be performed later on. The elementary volume of the first model of size (xmax~ Xmin» T) i n t h e xt-space (see figure 6.2a) will

be used. Let the particle distribution in the x-space at time tQ be denoted by U(x;t_). This means that

there is an amount of U(x;tQ) particles at point x.

Figure 7 is a perspective sketch of the xtU-space. The particle distribution at time tQ-T is denoted by

U(x;tQ-T). A relation between these particle distribu­

tions can be derived as follows.

Fig. 7. Particle distributions at t and t -T.

o o

A particle starting from point 1 has the conditional

probability F(X,;T) of arriving' at point A.

Consider one single point A of the elementary volume with coordinates xA,t0 (see figure 7 ) . The particle

amount at point A is U(x^;t0). The particle amount at

point 1 with coordinates xj,tQ-T is U(x^;t0-T). This

amount, however, can also be denoted by U(x^-X^;t0-T),

where

X,= x.- x, (6.1)

is the distance between x. and x,. So 1 A

U(xl 5to-T) = U(xA-Xl 5to-T) (6.2)

Not all particles starting from point 1 would arrive at point A, after a given time interval T. They have the conditional probability of F(X,;T) of travelling a dis­ tance X, in the given time interval T. There is an

amount of U(x.-X,;t -T) particles starting from point 1. Only the part U ^ - X ^ t -T).F(X1;T) will arrive at

point A. In other words, particles starting from point 1 contribute an amount of U(x.-X,;t -T).F(X,;T) par­ ticles to U(x.;t ) . Similarly, particles starting from

A O

point 2 (see figure 7) contribute an amount of U(xA-X2;t -T).F(X2;T) particles to U(xA;t ) . Taking

into account all points 1,2, which correspond to the set D(X;T) (5.1.1) of all possible displacements, the total amount of particles found at point A will be

The point A with coordinate x, at time t has been A o chosen arbitrarily. So the subscripts A and o can be omitted, giving a general relation for arbitrary values of x and t,

U(x;t) = E U(x-X;t-T).F(X;T) (6.4) D(X;T)

where X is an element of D(X;T).

Such a relation could not be found if we had used the elementary volume of figure 6.1. The elementary volume of figure 6.2b, however, is usable. This would lead to

U(x;t+T) = E U(x-X;t).F(X;T) (6.5) D(X;T)

which is equivalent to (6.A). This can be shown direct­ ly by replacing t by t-T in (6.5). We will proceed with

(6.4) since the unknown function U(x;t) is written ex­ plicitly in that equation.

If the particle distribution u(x;t) is continuous, then (6.4) becomes

u(x;t) = £ u(x-X;t-T).F(X;T) (6.6) D(X;T)

Expression (6.6) is not meaningless. For, a continuous particle distribution could be translated, for instance, over a discreet distance X during a time interval T. If X is also continuous, i.e. X is an element of a con­ tinuous set D(X;T), and the corresponding conditional

probability density function is denoted by f(X;T), then (6.6) becomes

u(x;t) = f u(x-X;t-T) f(X;T) dX (6.7)

D(X;T)

For the sake of simplicity, we only proceed with ex­ pression (6.4) or (6.7). These equations have been de­ veloped for an arbitrary small value of T. So they are still valid if T tends to zero.

U(x;t) = lim £ U(x-X;t-T) F(X;T) (6.8) T-fO D(X;T)

u(x;t) = lim ƒ u(x-X;t-T) f(X;T) dX

T . n J (6.9)

T^ ° D(X;T)

These equations denote the change of the particle dis­ tribution with respect to an elementary volume men­ tioned in section 5. So they can be considered as the governing equations of the main problem. Expression (6.7), in particular, is an integral equation with the unknown function u(x;t) and a kernel f(X;T). Such an integral equation is not that easy. We know nothing about the kernel f(X;T). Even if the kernel f(X;T) is given, it has to satisfy some conditions for the solv­ ability of the integral equation. Yet we will try to find a solution of (6.7) in very simple cases (see sec­ tion 7 ) .

The kernel f(X;T) is, in many cases, not given explicit­ ly. In such a case the integral equation cannot be

treated at all, at least not in its original form. An­ other approach will be required (see section 8 ) .

Similarly, by considering an elementary volume of size (T - T . , X) of the second model in the xt-space,

v max m m '

we will obtain the following integral equation

u(t;x) = ƒ u(t-Tjx-X) f(T;X) dT (6.10) D(T;X)

where the functions u(t;x), f(T;X), and the domain D(T;X) should be defined accordingly.

In the third model, we can say, as mentioned before, that

(6.11)

such that the integral equation degenerates into the trivial

u(x;t) = u(x-X ;t-T) (6.12)

which physically means a pure translation of the par­ ticle distribution function. As a matter of fact, this model is only interesting if the displacement depends on where or when the particle is released. In that case, however, we are dealing with deterministic problems in statics or dynamics, which have their own, more effec­ tive, solving methods.

In the fourth model, the step parameter plays the role

f(X;T) =• ( 1

(o

for X = X o for X 4 XQ

of time. This leads to

u(x;p) = f u(x-X;p-P) f(X;P) dX (6.13)

D(X;P)

or, with the time behaving like a place coordinate,

u(x,t;P) = JJ u(x-X,t-T;p-P) f(X,T;P) dX dT (6.14)

D(X,T;P)

By the same reasoning, the method can be extended to more dimensions. In n-dimensional space, in which a point is denoted by the vector £ and a particle dis­ placement by the vector .X, the governing integral

equation for the first model will be, similar to (6.7),

u(£;t) = f...fu(x-X;t-T) f(X.;T) dX. (6.15)

D(X;T)

where the n-dimensional domain D(X.;T) should be defined accordingly. The abbreviated notation dX_ stands for

dX,dX„ dX . This notation will be used frequently without further explanation.

The other models in more dimensions can-be treated similarly.

7. EXAMPLE SATISFYING THE INTEGRAL EQUATION

7,1. One-dimensional

Using the first model in the one-dimensional space, the problem is governed by the integral equation (6.7). The displacement probability field is supposed to be

Gaussian, and the conditional probability density func­ tion is given by

2

f ( X ; T ) = - 4

exp{- <

"ffl } (7.1.1)

with the parameters

y(X) = X = AT (7.1.2)

02(X) = 2BT (7.1.3)

In this context, A and B can be called the field con­ stants. By (7.1.3) it is obvious that B must be great­ er than zero. This is one of the cases in which the

'elementary volume' becomes infinite. The size in t-direction equals T and remains finite, but the dis­ placement X in x-direction during that finite time in­ terval T can vary from -co to + <», theoretically. This is physically impossible, but it is caused by the im-perfectness of the concept. It also happens in diffu­ sion problems. Nevertheless, it is accepted prac­ tically, because the probability for a particle of undergoing a large displacement tends to zero.

The constants A and B have a particular physical inter­ pretation, as we will see later on.

The governing integral equation becomes, according to ( 6 . 7 ) ,

CO

u(x;t) = f u ( x - X ; t - T ) f(X;T) dX (7.1.4)

_ 0 0

the kernel of which is denoted by ( 7 . 1 . 1 ) .

Now, what will be the particle distribution density function u(x;t) if we have a point injection at the origin x=0 ?

A point injection can be described as releasing simul­ taneously a sufficiently large number of particles at a given point. Since the amount of particles is normal­ ized, i.e. made equal to unity, we can formulate the initial condition as

u(x;0) = 6(x) (7.1.5)

in which the Dirac delta function 6(x) has the follow­ ing properties:

e

6(x) = 0 for x 4 0, and lira A ( x ) dx = 1 (7.1.6)

Its boundary condition in an unbounded space can be written as

lim u(x;t) = 0 at any time t (7.1.7) x-»-±<»

This example can be solved by means of Fourier Trans­ forms. Let the Fourier transform of u(x;t) with res­ pect to x be denoted by u(£;t). Then we can write (see, for instance, ref. 14, p. 19)

00

u(C;t) = /*u(x;t) exp(i£x) dx (7.1.8)

2

(with i = -1), and reversely,

u(x;t) = ^ T u U ; t ) exp(-ixe) d£ (7.1.9)

Similarly,

?(C;T) = ƒf(X;T) exp(i^X) dX (7.1.10)

The right hand side of (7.1.10) is known as the charac­ teristic function or the mathematical expectation

ê (i£X) of the probability density function f(X;T).

Substitution of (7.1.1) in (7.1.10) gives a well-known integral with result (see, for instance, ref. 12, p. 790),

f(C;T) = exp(iAT£-BT£2) (7.1.11)

u(£;t) = f /*u(x-X;t-T) f(X;T) exp(iCx) dX dx

-°° -°° (7.1.12)

According to the convolution theorem (see, for instance, ref. 14, p. 2 5 ) , it follows from (7.1.12):

u(5;t) = u U ; t - T ) F(?;T) (7.1.13)

Now (7.1.13) can be considered as a difference equa­ tion in the unknown function u with respect to the dis­ creet step T. This equation is subject to an initial condition which can be obtained by transforming the original one (7.1.5),

oo

u(5;0) = f ó ( x ) exp(i£x) dx = 1 (7.1.14)

—oo

For the moment, we cannot yet substitute t=0 in equa­ tion (7.1.13), but we can start with t=T, getting

u(£;T) = u(S;0) r(e;T) (7.1.15)

Taking into account (7.1.11) and (7.1.14), expression (7.1.15) becomes

u U ; T ) = f(C;T) = exp{(iA^-BC2)T} (7.1.16)

If we now substitute t=2T in equation (7.1.13), we ob­ tain

u(£;2T) = u(S;T) f(5;T) = f2(£;T) = f(S;2T) (7.1.17)

Apparently, it applies for any number n:

u(S;nT) = fn(5;T) = f"(£;nT) (7.1.18)

as the exponent in (7.1.16) is linear in T.

Now the step T can be taken arbitrarily small. The step T can be such chosen that any required t can be made equal to nT,

t a nT (7.1.19)

Combining (7.1.18) and (7.1.19) gives the solution of equation (7.1.13):

GU;t) = HE.it) (7.1.20)

Then, from (7.1.20) and (7.1.1), it follows directly

» ( x ; t ) = - b e x p { - . Ü ^ } (7.1.21)

We come back to a previous remark that it was not pos­ sible to substitute t=0 directly in the difference equation (7.1.13). Taking t=T with T in the limit tend­ ing to zero, will resolve this difficulty.

It can be checked easily, that expression (7.1.21) satisfies the initial condition (7.1.5) and the bound­ ary condition (7.1.7). It can also be shown that (7.1.21) fulfils the integral equation (7.1.4). The

latter will be performed as follows. Substitution of (7.1.1) and (7.1.21) into the right hand side of (7.1.4) gives ^o 2AB(t-T) { WKZ l ) }

*ik w s s * }

x / L / (x-X-A(t-T))

_ (X-AT)

1

T h e integral (7.1.22) c a n be worked o u t . On substitut­ ing

X - A T = q + dX = dq (7.1.23)

the integrand o f (7.1.22) becomes

e x p { - ( a q2- 2 b q + c ) } (7.1.24)

in which

a = - • b o x"A t • c - (x"At) (7 1 25}

The expressions a,b,c of (7.1.25) are similar to those of de Josselin de Jong (ref. 9, p. 24), arrived at, however, by a different method. Now (7.1.24) can be decomposed as

2

exp {-{(/a.q - ^-)2+ (c - £ )}} (7.1.26)

and after substituting

/ b , dv 2/BT(t-T) , ,-. , „-,,. / a , q" ~ ^ '= v * d q = 7a" = A (7.1.27) t h e i n t e g r a l ( 7 . 1 . 2 2 ) becomes 2 °° u ( x ; t ) = — - — exp(-c + - ) / * e x p ( - v2) dv ( 7 . 1 . 2 8 ) 2WÏÏT a Ji

Since (7.1.28) contains a well-known integral

oo Z"exp(-v2) dv = /IT (7.1.29) — 0 0 it becomes 1 h2 u(x;t) = — - — exp(-c + - ) (7.1.30)

2/rrËF

The term between brackets of (7.1.30) can be developed by (7.1.25), giving the expected result

2 2

_c . b _ (x-At) ( 7 j 3 n

c + a - A B t W . I . J I ;

This forms the last link of the proof, as (7.1.30) be­ comes identical with (7.1.21).

It has been mentioned before, that the use of the second model consists of interchanging the roles of the place coordinate and the time variable and at the same time the displacement and the time interval. Apart from its physical interpretation, the procedure

is exactly the same.

Obviously, the third model cannot be used here. The fourth model leads, similarly, to the solution

2

u ( x ; p ) = — ^ exp {- (*:$ } (7.1.32)

2/HEp

'

satisfying the integral equation (6.13). This solution, however, is not usable, as it does not contain the time variable t. For, even in the one-dimensional case, one requires at least two variables, place x and time t. It might be usable, if the step parameter is, for in­ stance, proportional to the time variable t. But then (7.1.32) is, apart from a constant proportionality factor, exactly the same as (7.1.21). Therefore, it would be better to treat this in the next subsection.

7.2. Two-dimensional

Let in two-dimensional space the displacement probab­ ility field be Gaussian again. Its parameters are the

average values y(X) and y(Y), the variances O (X) and

2

O (Y), and the covariance cov(X,Y). The correlation

coefficient can be denoted by p(X,Y), which is defined

as cov(X,Y)/a(X)0(Y). Usually, the probability density

function is written as

f(X.Y) =

1

2TTG(X)O(Y

)/!-(?

exp

* 2(l-p

)

jiXzMmii _

^-KX)HY-,(Y)}

IMQiijI

<T(X) °

cT(Y) '

in which p stands for p(X,Y).

For our purpose, however, it will be more convenient

to write (7.2.1) in another way (using the first model),

f(X,Y;T) -

1

4UT/|BT

, b

(X-a

T)

-2b

(X-a

T)(Y-a

T)+b

(Y-a

T)

>

x exp

^ ^ J

4 |B| T

(7.2.2)

where the vector |i and the symmetric matrix B, with

the given time interval T, are the new parameters.

They are defined as

a =

; B =

b

ll

12

12

22

(7.2.3)

IBI = det B = b

,b

- b

In (7.2.2) |B| must be greater than zero. The case |B| < 0 falls outside the scope of this study, and the case |B| = 0 will be discussed in section 13. Since (7.2.1) and (7.2.2) are identical, there exist the following relations between the parameters, which can be verified easily.

' U(X) . U(Y) . = ' X ' . Y . = al ■ a2 ■ T = a T (7.2.4) ' a2(X) cov(X.Y) cov(X.Y) 02(Y) f bll b12 b12 b22 2T = 2BT (7.2.5) pf x V) = cov(X,Y) MU.ï; - ö(X)o(Y) '12 /bllb22 (7.2.6)

The governing integral equation in an unbounded space is, similar to (6.15),

00 00

u(x,y;t) = /7u(x-X,y-Y;t-T) f(X,Y;T) dX dY (7.2.7)

the kernel of which is given by (7.2.2). Now consider a point injection at the origin (x=0, y=0). Its initial and'boundary conditions can be formulated as

lim u(x,y;t) = O at any time t and for any y

x-y±<*>

lim u(x,y;t) = 0 at any time t and for any x

y-V+oo

Similar to (7.1.6) the delta function <5(x,y) has the following meaning:

<5(x,y) = 0 for x^O or y?ÉO

e (7.2.9)

lim //S(x,y) dx dy = 1

Making use of a double Fourier transform, similar to the method used in subsection 7.1, the solution of

(7.2.7) will be

u(x

'

4ÏÏE7JBT

x exp

{-2 {-2 b2 2(x-a1t) -2b12(x-a1t)(y-a2t)+bn(y-a2t) x

4 |B| t > (7.2.10)

This already satisfies the initial and boundary con­ ditions (7.2.8). We still have to show whether (7.2.10) fulfils the integral equation (7.2.7). It does. The working out is not written here, because it is elaborate. Essentially, it does not differ much from 7.1. Moreover, the more-dimensional case will be treated completely.

to the use of the second and third model.

The method used in this two-dimensional case, denoted by the place coordinates x and y and the time variable as parameter, can be applied for the one-dimensional fourth model, in which the place coordinate x, the time variable t and the step p as parameter, play the leading part. Let the field be Gaussian with parameters U(X), u(T), 02( X ) , a2(T), p(X.T). The correlation co­

efficient p(X,T) is defined as cov(X,T)/{o(X)a(T)}. Then, similar to (7.2.1), the probability density function reads

f(x.T) * . . exp{-2nè^y{

2ïïö(X)0(T)/ï^pT l A i P J

<{X-u(X)}

{X-u(X)}{T-u(T)} {T-u(T)}

,)

I

2

" * > 0(X)0(T)

2

1/

(7.2.11) where p stands for p(X,T). For our purpose, similar to (7.2.2), expression (7.2.11) can be written as

f

<

»

'

> - W l B l

, b2 2(X-a1P)Z-2b1 2(X-a1P)(T-a2P)+bn(T-a2P)2 ,

X eXP

\ 4lBl~P /

(7.2.12)

where the vector a_ (with components a, and a„) and the symmetric two by two matrix B (with elements b,,, b,„= b2i, and b„2) and the given step magnitude P, are the

new parameters. Again, |B| = ^ii^^o" ^IO m u s t ^e greater than zero. Since (7.2.11) and (7.2.12) are identical, there exist the following relations

between the parameters, which can be verified easily.

' y(X) '

.

u(T)

.

It should be noted that (7.2.13), (7.2.14), (7.2.15) are similar to (7.2.4), (7.2.5), (7.2.6), respectively, but they are not the same. The governing integral

equation in the xt-space, using the fourth model, is (6.14) with (7.2.12) as kernel.

Now consider a point injection at the origin (x=0, t=0) of the xt-space. The time variable will be treat­ ed as a place coordinate, while the step parameter p now plays the role of time. The 'initial' condition then can be written as

u(x,t;0) = 6(x,t) (7.2.16)

where the delta function has the same meaning as in (7.2.9). The boundary conditions read

lim u(x,t;p) = O for any t and any p

X > ±~ (7.2.17)

lim u(x,t;p) = 0 for any x and any p t-*-±°°

It should be remarked that the physical interpretation of these conditions cannot be imagined easily. An a t ­ tempt is perhaps by considering the step parameter as something like a dimensionless progress in space or time. Now the solution, similar to (7.2.10), will be

u ( x

'

4VIBI

r

AIBTP

/

This 'solution', however, still contains the parameter p . The final solution may not contain p anymore. This subject has been studied before by de Josselin de Jong (ref. 9 ) . A large part of this subsection has been taken over from his study. The required particle dis­ tribution function in the xt-space can be obtained by integrating (7.2.18) over the p-domain, i.e. over all possible values of p . So we get

00

u ( x > t ) =

4ÏÏ71ËT ƒ

"

" p

^ (7.2.19)

0

1 2 2

a = 4 IBI' ^b22al ~ 2 b12ala2+ blla2^ (7.2.20)

B = 2 l B l{ ( b2 2ar b12a2) x + ( bl la2 - bi2 al) t } C7.2.21)

Y = JJQ\ (b22x 2 _ 2 b12x t + bl lt 2 ) (7.2.22)

The integral (7.2.19) can be worked out as follows (ref. 1, 9, and 22). The f i r s t step i s taking exp(g) outside the integral sign, as $ does not contain p.

oo

u ( x > t ) = A^TfBl" e x p ( 3 ) /e xP < - aP " p> ^ (7.2.23)

2 2 By substitution of ap = y and cty = Q t n e integral of (7.2.23) becomes

ƒexp(-ap - J ) ^ = 2 ƒ e x p ( - y2- \ ) ^ = 2 I (7.2.24) o o y

Differentiation of I from (7.2.24) with respect to q yields 2 g = - 2 q fexp^y2- A-) ^ (7.2.25)

o

Differentiate K from (7.2.26) with respect to q,

2

jj| = - 2q ƒ exp(-z2- % ) ^f = -2q I (7.2.27)

From (7.2.26) and (7.2.27), it follows

2

^-4 + ^4^-41 = 0 (7.2.28)

dq2 q d q

By substitution of q = y r equation (7.2.28) becomes 2

r2 4 - | + r il _ r 2ï = o (7.2.29)

dr

The solution of equation (7.2.29) can be written di­ rectly,

I = Cj_ Io(r) + c2 KQ(r) (7.2.30)

in which I and K are modified Bessel functions of o o

order zero, while c, and c„ are arbitrary constants to be determined by boundary conditions. From the con­ dition that I = 0 if r goes to infinity, it follows c, = 0, and (7.2.30) becomes

I = c2 KQ(r) (7.2.31)

For small arguments, (7.2.31) can be approximated by (see, for instance, ref.1, p.375)

I - c

K

(r) - i c

(7.2.32)

On the other hand, it follows from (7.2.24),

OO ry

I =ƒexp(-y

-

2) ^ (7.2.33)

0

which for small values of q (or r) can be approximated

by

OO OO

I = ƒ exp(-y

) &.

I y

p ( - X ) X

dX =

\ T(0) = \

0 0

(7.2.34)

in which

T denotes the gamma function. Now, it follows

from (7.2.32) and (7.2.34), c„= 1, and hence

I = K

(r) = K

(2q) =

(7.2.35)

F i n a l l y , ( 7 . 2 . 2 3 ) becomes

" (x. t ) = ^ 7 T 5 T exp(3) K J 2 v ^ ) ( 7 . 2 . 3 6 )

2TT/|B _o

where K denotes the modified Bessel function of the

o

second kind and order zero. If x and t are sufficient­

ly large, then v^xy ^ 1» and the Bessel function can

be approximated by (ref. 1, p. 378)

C (2^0Y) = \/—^- exp(-2v/oY) (7.2.37)

By substituting (7.2.37) into (7.2.36) the particle density in the xt-space becomes

u(x,t) = * exp(B - 2v^OY) (7.2.38) 4yTr|B|v^Y

The particle density u(x) at a particular time t can be obtained (ref. 12, p.770) by substituting t = t into (7.2.38),

u ( x )t = t = X u(x,to) (7.2.39)

o

where the factor X must be chosen such that

CO

X fu(x,tQ) dx = 1 (7.2.40)

— 0 0

7.3. More-dimensional

In n-dimensional space, it will be more convenient to use matrix and vector notation. Let the vector J_ be

the stochastic displacement of a particle in that space during a given time interval T. And suppose the displacement probability field to be Gaussian again. Then the conditional probability density function can be written as

i / (X-aTVB ^ X - a T ) .

f(X;T) = exp { J

(2/ÜT)n/|ï]" { 4 T '

in which the n-dimensional vector a_ and the n by n sym­ metrical matrix B , with the given time interval T , are the given parameters. The accent means transposition. A restriction is again: |B| = det B > 0.

Similar to (7.2.4) and (7.2.5) the average displacement X_ and the covariance matrix C of the displacement are, respectively,

X = a. T C = 2 B T (7.3.2)

Let x. denote an arbitrary point in this space, then the governing integral equation can be expressed by

OO 00

u(x_;t) = [...[ u(jc-X;t-T) f(X;T) dX (7.3.3)

_0O —CO

n-fold

with (7.3.1) as kernel. A point injection at the origin >c = 0_ of an unbounded n-dimensional space can be de­

scribed by the initial condition

u(x.;0) = 6(x) (7.3.4)

with

e e

6(x) - 0 for x 4 0 and lim f.. ./*6(x) dx. = 1 (7.3.5)

e

+°-e

-{

n-fold and the boundary condition

lim u(x.;t) = 0 at any time t (7.3.6) x-+-±°°

Now a multiple Fourier transform with respect to x. and X_ can be applied (see, for instance, ref. 14, p. 43),

00 00

u(£;t) = C...C u(x.;t) exp(il'x) dx (7.3.7)

— 0 0 _ 0 O

n-fold

where the accent denotes a scalar multiplication of the two vectors £ and x.« Reversely,

00 OO

u ( x j t ) - — L - f...fu&t) exp(-ix'i) di (7.3.8)

(2TT)" y -/

— o o —oo

n-fold

By (7.3.7) the integral equation (7.3.3) becomes

oo oo

u(£;t) = C...C u(x.-X_;t-T) f(X;T) exp(i£'3c) dX dx

2n-fold (7.3.9) and, according to the convolution theorem (see, for in­

stance, ref. 14, p. 2 5 ) , expression (7.3.9) reduces it­ self into the following difference equation,

u(£;t) = u(£;t-T) ?(£;T) (7.3.10)

where f denotes the Fourier transform of (7.3.1) with respect to X_. The unknown function u(£;t) has to sa­ tisfy the initial condition

u(£;0) = 1 (7.3.11)

equals one. By a similar reasoning as in the one-dimen­ sional case, see subsection 7.1, the solution of

(7.3.10) reads:

u(l;t) = f"(l;t) (7.3.12)

Then, with (7.3.1), we get

x i (x-ajO'B'^x-atV

(2/f?t)n/fBT { 4 t '

(7.3.13) The solution (7.3.13) satisfies the initial condition

(7.3.4) and the boundary condition (7.3.6). It will be shown that (7.3.13) also satisfies the integral equa­ tion (7.3.3). Substitution gives

OO 00

u(x/,t) = _ ƒ. . . ƒ exp <

{4uA(t=T)}n|B| J J {

1 1 _ 0 0 —OO

n-fold

i {xrX-a.(t-T)}'B"1{xrXra.(t-T)} (X.-aTyB^QC-aT),

{ 4Ü=T) — « ; d*

(7.3.14) By introducing

£ = X - aT -»• d£ = dX (7.3.15)

the integrand can be written as

+ 2(x-at)'B Xa 4 ( tj ;T ) - (x-at)*B_1(x-£t) A(t=T)}

(7.3.16) Let, for shortening,

k = 4T(t-T) a n d — = 4(t-T) ^—~—t^ (7.3.17)

then the integrand becomes

exp {- ka/B"1^ + 2 s.'B-1£ - 4(t-T) s/B_1s. I (7.3.18)

which can be decomposed into

s s

exp {- ( A CL- ^ ) ' B_ 1( A 2r ^ ) + f £ -4(t-T) } s'B_1s }

(7.3.19) And the next substitution

s^

A £ - 7£ = v. ->• d£ = (A)~ndv. (7.3.20)

brings the integrand into the form

exp|- x'B"1!- ^ " ^ ^'B"1^! (7.3.21)

Then expression (7.3.14) becomes

u ( x ; t )

1 —

(_ Mt£üi

.

-i

\

00 CO

X ƒ*.. ƒ exp(-v_'B *v) dv. (7.3.22)

f— oo —00

'n-fold

By (7.3.17) the exponent outside the integral sign turns out to be

(x-at)'B (x-at)

exp

_{{" 4t /}

_'

_'

The last step is the computation of the integral from (7.3.22). In the n-dimensional space, v/B v_ = constant

represents an n-dimensional ellipsoid. For, det (B ) is greater than zero, which follows directly from the assumption that det B > 0. The coordinate axes v^ can be rotated, by an orthogonal rotation matrix R, such that the principal axes of the ellipsoid fall along the new coordinate axes v .

v = R v_ (7.3.24)

The integrand becomes

exp(-v/R*B-1Rv_) = exp(-v/Dv.) (7.3.25)

in which D = R'B R has become a diagonal matrix with positive diagonal elements. Actually, this part is an eigenvalue problem. The corresponding transformation is called a principal axes transformation (ref. 12, p. 62) or a unitary transformation (ref. 5,p. 18-4). The

an orthogonal matrix. This has a physical interpreta­ tion, namely that by this rotation the volume of the ellipsoid does not change. The integral becomes

CO U U /• • •/ exp(-VDy_) dv_ = _ o o - c o n - f o l d oo oo i = n /.../expj- £ dü vi } d vid v2 d vn (7.3.26) — C O — C O 1 — 1 n-fold

Now the n-dimensional ellipsoid can be transformed into a unit sphere, by introducing new coordinates v* which coincide with the old ones v^ » but possessing other scale factors.

v, = - Z Z w. (i « 1,2 n) (7.3.27)

/ di i

The Jacobian | j | of this transformation then equals

i=n

i=

i T d T /]DT •|R«B-

R| /fFH"

(7.3.28)

and the integral becomes

i=n

/ • /

e x

p { ^ E "I} "W

d

V

w

2

d

w

n

= /lBTi

n

_ o o —co 1 = 1

n - f o l d

The integral I should be computed separately. This can be performed easily, since I can be written as a pro­ duct of n single integrals.

i=n /. „

I - P I /

P(-

-i) dw. = ( A ) " (7.3.30)

i=l J_ 1

Now the proof can be completed by substitution of these results in (7.3.22), yielding

i i ( x - a t y B ' ^ X r a t ) * ^

u ( x ; t ) = ±- exp { 7- f / ^ T ( A )

(2Tr/t)

|B|

'

(7.3.31)

which is exactly the same as (7.3.13).

This was the use of the first model. The same remarks as in subsection 7.2 apply to the second and third model. The use of the fourth model in n-dimensional

space leads to the integral equation, similar to (6.13) and (6.14),

u(x.;p) = ƒ••/u(x-X;p-P) f(X;P) dX. (7.3.32)

D(X;P)

It should be noted that £ here is an (n+l)-dimensional vector with the place coordinates as the first n compo­ nents and the time variable t as the last component. A similar remark applies to the vector X_. Let the field be Gaussian,

. . (X-aP)'B ^ X - a P ) ,

f(I;P) = V = = e x p { - TO ) (7.3.33)

(2/SP)

/|¥T

'

in which the (n+l)-dimensional vector and the (n+1) by (n+1) symmetrical matrix B, with the given step magni­ tude P, form the parameters, similar to (7.2.13) and (7.2.IA),

I = a _ P C = 2 B P (7.3.34)

If we have a point injection in the (n+l)-dimensional xt-space, the initial and boundary conditions of which are the (n+l)-dimensional equivalent of (7.2.16) and (7.2.17), then the solution of (7.3.32) will be, similar to (7.2.18) and (7.3.13),

I i (£-ap)'B~ (£-£P)\

u(x.;p) = exp \- T- ) (7.3.35)

(2/^)

/ÏBT *

P '

Now ( 7 . 3 . 3 5 ) should be i n t e g r a t e d over t h e p-domain.

! r i ( x - ^ ' B - ^ x - a p ) u(x) = ƒ exp < - ; > —K— " (2 V/?)VTBT { < 4 p ' (/P)n (7.3.36) By i n t r o d u c i n g a, 3 , y, s i m i l a r t o ( 7 . 2 . 2 0 ) , ( 7 . 2 . 2 1 ) , and ( 7 . 2 . 2 2 ) , ( s e e a l s o r e f . 9, p . 2 4 ) , a = \a!?>~la_ ( 7 . 3 . 3 7 )

3 = ia'B_ 1x (7.3.38)

Y = i x ' B_ 1x (7.3.39)

the integral (7.3.36) becomes

uu

u(x) - ~- fexv(-ap + B - J) - ^ — (7.3.40)

(2/ir)n/ÏBT(/ P (/P)n

The factor exp(3) can be taken outside the integral sign, as it does not contain p. The remaining integral,

oo

J = /"expC-ap - h - & - (7.3.41)

0J P (/P)n

can be evaluated by a similar method used in subsection 7.2. It is not treated here, because it is elaborate. Besides, it is not essential in this study. For the moment, we can conclude cautiously that it does not pay to use the fourth model if the integration over all the possible steps cannot be executed easily.

8. DERIVATION TO DIFFERENTIAL EQUATION

Using the first model in one-dimensional space, with continuous displacements, the problem is governed by the integral equation (6.7). The kernel f(X;T) of this equation, however, is not given explicitly. So the in­ tegral equation cannot be treated directly, and we have to look for another method. The integral equation has been derived for an arbitrary small value of T. So it is still valid if T tends to zero.

u(x;t) = lim f u(x-X;t-T) f(X;T) dX (8.1)

l ^ D(X;T)

If T goes to zero, the displacement X goes also to zero (5.1.4). Then the elementary volume becomes infinite-simally small, as always required in developing a dif­ ferential equation. Now it is reasonable to assume that the function u(x-X;t-T) can be developed in a Taylor series,

u(x-X;t-T) = u - X ux- Tut+ |, (X2ux x+2XTux t+T2ut t)- ....

(8.2)

where u stands for u(x;t) and a variable as subscript means partial differentiation with respect to that variable. On substituting (8.2) in the integral equa­ tion (8.1), and integrating term by term, we are con­ cerned with integrals of the form

These integrals are the moments u or the expected values € of X (i = 0,1,2, ). In this study they

can be identified with average values X1 (denoted by an overbar).

U±(X) - ëtX1) = X1 = ƒ X1 f(X;T) dX (8.A) D(X;T)

So, when the function f(X;T) is not given explicitly, we wijl do with its moments (8.A), if they exist. The moment of order zero equals one, according to (5.1.3). Now (8.1) becomes

u = limiu-(Xu) -Tu + 4.{(X2u) +2(XTu) „ + T V J - ..A

m n I x t 2'. 'xx 'xt tt )

(8.5) Since u cancels in (8.5), dividing by T gives

If the displacement X tends to zero, all its moments X (i = 1,2, ) also go to zero. Then (8.6) becomes

— 2 3

V £jj {-(l u )x+ Ï A U)xx" l!(f U)xxx+ - • } ( 8'?)

It has been assumed that the displacement X does not depend on where (x) or when (t) the particle has been released. So (8.7) becomes

/ X 1 X^ 1 X^ ) u = lim < - ~ u + 77, i: u - T . ™ u + .. •> (8.8) fc r^o \ T x 2! T xx 3!T xxx f or u = - A u + B u - C u +... (8.9) t X XX XXX

with the coefficients A,B,C,... defined as

— 2 3 X I X I X

A = lim =r ; B = lim ■=-, ™ ; C = lim TT, „■ ; etc. (8.10) T O T-K) T O

If the limit values (8.10) exist, then (8.9) is mean­ ingful, and we can say cautiously that the particle spreading process is governed by the partial differen­ tia! equation (8.9). The coefficients A,B,C,... have a physical meaning. A is the average displacement per unit time (see (8.10)). So it can be interpreted as the average velocity of the particles. B is a measure of the deviation of the particles from the average value. The higher order moments are rather difficult to interpret. At any rate, they all define the con­ ditional probability density function f(X;T) implicit­ ly. Since T tends to zero, the n-th order (n 2 2) moments can be approximated by the reduced moments

(X-X)n

. This can be shown as follows. Newton's binomial theorem gives

i=n n-i

(X-X)n = £ (J) X1 X (8.11)

We have seen from (5.1.3) that X°= 1 and from (8.10) that X1 (i = 1,2,...) should be of order T. Then, only the last term Xn of (8.11) is of order T, while other terms are of order T^ or higher. So (8.10) can be re­ placed by

X (X-X)^ (X—X) A = lim TZ ; B = lim 9I T' ; C = lim \,.r/ ; etc.

T-+0 T-K) T-K) (8.12) Equation (8.9) is of infinite order. Such an equation is meaningless, as well mathematically as physically. In real problems there should be always physical

restrictions such that the equation can be reduced to, or at least approximated by, a finite, practically low, order equation (see section 9 ) .

The use of the second model leads to the equation

V - K ut + L ut t- M ut t t + ... (8.13)

in which

T 1 T^ 1 T^

K = lim rr ; L = lim r. - ; M = lim -r. TT- ; etc. (8.14)

X-K) A X-K) Z- A X-K) J' A

The third model cannot be used properly. If one would use this model at any cost, then it would lead to the equations,

instead of (8.9) or (8.13). The solution of (8.15) is, respectively,

u(x;t) = u(x-At) and u(t;x) = u(t-Kx) (8.16)

They denote physically a pure translation of the

particle distribution function along the x- or t-axis. The use of the fourth model in one-dimensional space, similar to (8.9) and (8.10), leads to the equation

u = - A u + B u - C u + . . . (8.17) p x xx xxx in which X I X 2 1 X"^ A = lim =■ ; B = l i m ■=• tj ; C = lim ■=-. — ; e t c . ( 8 . 1 8 ) P-K) r P->0 Z- V P-K) J' V

The described method can be extended to more complicat­ ed cases. Place- and/or time-dependent fields fall out­ side the scope of this study. An extension to more-dimensional cases will be performed now.

In two-dimensional space the integral equation will be

u(x,y;t) = lim f f u(x-X,y-Y;t-T) f(X,Y;T) dX dY

D(X,Y;T) ( 8 > 1 9 )

On developing u(x-X,y-Y;t-T) in a Taylor series,

u(x-X,y-Y;t-T) = u - Xu - Yu - Tu +

X y t (8.20

and integrating term by term, we obtain

ut- - a l V a2uy + bnux x + 2 b1 2ux y + b2 2uy y- ... (8.21)

where

X Y a,= lim ™T ; a„= lim „

1 T-K) T-K)

,. 1 Xz , . 1 (X-XV

b,,= lim 9"t ^r = lim ö"i T

1 1 T-K) T-K) r 1 XÏ b10= lim -x. -sr 1 Z T-K) lim -y, T-K)

I IX:

:-XHY-Y)

By considering ^ « ^ o a s t'ie components of a vector a_,

am B,

and b,,,b..„,b„2 a s t n e elements of a symmetrical matrix,

a = B =

bn b1 2

b12 b22

(8.23)

equation (8.21) can be written as