Marco Scianna and Matteo Icardi

An Introduction to Mathematical Modeling Relazione per il corso di Equazioni della fisica matematica

Gennaio 2006

Contents 1 Mathematical Models and Modeling Scales

5

1.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2

Definition of Mathematical Model

. . . . . . . . . . . . . . .

6

1.3

Modeling Scales and Representation . . . . . . . . . . . . . .

8

1.3.1

Passages from one scale to the others . . . . . . . . . .

11

1.4

Elementary Examples of Models at Different Scales . . . . . .

14

1.5

Critical Analysis . . . . . . . . . . . . . . . . . . . . . . . . .

17

1.5.1

Boundary Conditions

. . . . . . . . . . . . . . . . . .

17

1.5.2

Complexity . . . . . . . . . . . . . . . . . . . . . . . .

21

2 Models at Macroscopic Scale

23

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

2.2

Deduction of Models . . . . . . . . . . . . . . . . . . . . . . .

24

2.3

General deduction and different type of flux . . . . . . . . . .

27

2.3.1

Vehicular Traffic Model . . . . . . . . . . . . . . . . .

30

Classification and Formulation and diffusion model properties

32

2.4

3 Discretization Methods of Macroscopic Models

35

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

3.2

Discretization of Diffusion Models

. . . . . . . . . . . . . . .

35

3.2.1

The finite difference methods. . . . . . . . . . . . . . .

36

3.2.2

The collocation method. . . . . . . . . . . . . . . . . .

39

Simulations of Diffusion Models . . . . . . . . . . . . . . . . .

46

3.3

4 Conclusions

51

3

Chapter 1

Mathematical Models and Modeling Scales 1.1

Introduction

Every physical system can be observed and studied to reach a knowledge of the inner structure and of its behavior. So traditional experiments are now supported by the construction of mathematical models of the physical system which have to reduce the complexity of the real system and to help its interpretation and prediction. Nowadays mathematical modeling plays a crucial role not only in physics but in all applied sciences (from natural to technological ones). This chapter is a general introduction of mathematical modeling science in which we give some definitions and tools. In particular: • Section 1.2 deals with the definition of mathematical model and of the mathematical quantities related to. • In Section 1.3.1 we report the definition of the different observation scales and the related representation methods and we give some ways to pass from one scale to the others. In this section there are also some other types of classification of mathematical models. • Section 1.4 shows some simple examples of mathematical models which correspond to the above classification. 5

6

1. Mathematical Models and Modeling Scales • In the last section there is a critical analysis of mathematical modeling sciences. We analyze the mathematical features of the model and we characterize a mathematical problem, also by studying its initialboundary conditions.

1.2

Definition of Mathematical Model

Let us present detailed definition of mathematical model. In order to do this we need some ingredients such as the following definitions: • Independent variables A variable is independent if it exists apart from the existence of the system. Usually independent variables are represented by time t ∈ [0, T ] ⊆ <+ and space x = {x, y, z} ∈ D ⊆ <3 , where [0, T ] and D are respectly the time interval and the volume in which we observe the system. • State variables The state variables is the finite dimensional vector variable u = u(t, x) : [0, T ] × D →
1.2 Definition of Mathematical Model

A Mathematical Model of a real physical system is an evolution equation, or a set of equations,such as u = u(t, x) : [0, T ] × D →
+ div(ψv + f ) = s,

where ψ is the quantity studied, ψv + f is the convective flux and s is the source term. If the variation is zero, we obtain a conservation equation. A Mathematical Problem is the mathematical model linked to its boundary or initial condition (for more details see Section 1.5) At this point we can introduce some important concepts. • Consistency A mathematical model is said to be consistent if the number of unknown dependent variables is equal to the number of independent equations. • Good formulation A mathematical model is said to be well-formulated if the evolution equation is associated with the correct number of initial or boundary conditions for its solutions.

7

8

1. Mathematical Models and Modeling Scales

• Well-posedness A mathematical model is said to be well-posed if it is characterized by existence, uniqueness and continuous dependence of the solution on the initial data. It can be proved that: wellposed =⇒ wellf ormulated but not the viceversa. A problem can be well-posed locally in time (∀t ≥ 0) or globally (∀t ≥ 0 and t −→ ∞). • Validity A mathematical model is said to be valid the information obtained from the model is sufficiently near to that obtained by experiments on the real system. In mathematic terms, introducing a norm for u , we can say that our model is valid if d=

ku-vk <ε kvk

where u is the state variable calculated with our modelization and v is a variable that expresses the real state of the system. There could be different choices of the norm (based on some kinds of average or an uniform norm) related to the needed approximation and it influences the above definition of distance.

1.3

Modeling Scales and Representation

Both observation and simulation of system of real world need the definition of suitable observation and modeling scales so that models and descriptions may correspond to different scales. The same system can be observed and described by the point of view of each singular particle or by using a larger scale considering suitable averages of the mechanical quantities linked to a larger number of particle. By the way the definition of small or large scale has a meaning which has to be related to the size of the object and ,to the volume containing them and to the homogeneous characteristics of the

1.3 Modeling Scales and Representation

system. In order to do this , let S be a real system, contained in a volume Ω set in a system of coordinate axes. We can suppose S as a set of elements S Ei , i = 1, . . . , n, i.e. S = ni=0 Ei . So we can give these definition: • Microscopic scale: a real system can be observed, measured and modeled at the microscopic scale if all single objects composing the system are individually considered, so that we make use of a mathematical model for each element Ei . • Macroscopic scale: a real system can be observed, measured and modeled at the macroscopic scale if suitable averaged quantities related to the objects composing the system are considered so we make use of a single model for the local average of the state of the elements Ei . The components of the state variable related to the macroscopic scale are called macroscopic observables. • Mesoscopic scale: a real system can be observed, measured and modeled at the mesoscopic (kinetic) scale if it is composed by a larger number of interacting individuals and the macroscopic observables can be recovered from weighted moments of the distribution function of the state of the system. Each of those types of representation lead to a different type of equation: • microscopic model is a global description of a system,so the variable state depends only on time and at this scale, the evolution equation is generally represented by ordinary differential equations (ODE) , like du = f (t, u; r). dt A mathematical model whose state variable does not depend on the space variables is also called discrete (otherwise it is continuous); • at the macroscopic scale u depends also on space because it is a local representation. Hence macroscopic mathematical models are generally stated in terms of partial differential equations, such as µ ¶ ∂u ∂ ku k = 1, 2, . . . ∂t = f t, x, u, ∂xk

9

10

1. Mathematical Models and Modeling Scales

There are also many other model classifications, the most important are the following. Classification by the state variable A mathematical model is dynamic if the state variable u depends on the time variable t, otherwise it is static. A system characterized by the reaction at an external force (ψ(t, u)orφ(u)) is called open else it is called close. Classification by the type of the state equations A mathematical model is algebraic if the state equation is an algebraic equation over the state variable and the set of parameters r characterizing the system f (u; t, x, r) = 0 where t and x behave as parameters in the same way as r. Using these kind of classification we can summing up with the following table: Static

Dynamic

Discrete

Algebraic equations

ODE

Continuous

PDE with respect to x

PDE with respect to x and t

Classification by the structure of the state equations A mathematical model is linear if the state equation can be written in the form Lu = 0 where L is a linear operator that satisfies the following condition L(u1 + u2 ) = Lu1 + Lu2 A mathematical model is nonlinear if the state equation can be written as Nu = 0

1.3 Modeling Scales and Representation

where N is a nonlinear operator that does not satisfy the former condition. A mathematical model is semi-linear if the state equation can be written as Lu + ²N² u = 0 where ² is a parameter, and the operator N² may depend on ².

1.3.1

Passages from one scale to the others

To show the ways to pass from one scale of representation to another let us give an example. Let us consider the motion of a fluid in a duct: at the microscopic scale, every particle composing the fluid is such as an element Ei of the real system S,so the motion of each particle can be described by Newtonian mechanics laws: ui = {xi , vi },

i = 1 . . . n,

where xi and vi represent position and velocity of the i-th particle. The passage from the microscopic to the mesoscopic scale is made introducing a statistical representation for the state of the system, making use of a statistical distribution f = f (t, x, v) ≥ 0 such that f dv = dn(t, x), if f is integrable, n(t, x) =

R R3

f (t, x, v) dv,

where n(t, x) is called numerical density and represents the number of particles for volume unity. Instead NΩ (t) =

R

Ω n(t, x) dx,

is the total number of particles in Ω at time t. If Ω is closed, this quantity to be constant: N (t) = N (t = 0) = N0 and we can normalize f to obtain:

11

12

1. Mathematical Models and Modeling Scales f (n) =

f N0

so that R Ω×R3

f (n) dxdv = 1.

f (n) dxdv = P(t, x, v ∈ [(x, v), (x + dx, v + dv)]) , is called probability density and {x, v} is the phases space. We can observe that f (n) is a probability distribution of microscopic states. Now let us introduce • the density mass, which is: ρ(t, x) = m n(t, x), where m represents the mass of the system. • and the flux, defined by: q(t, x) =

R R3

vf (t, x, v) dv,

or equivalently q(t, x) = n(t, x)U(t, x), where U = E[v] =

R vf (t,x,v) dv RR3 3 R f (t,x,v) dv

is the mean velocity. If v represents

the microscopic velocity, U is the macroscopic one. We now define the mean energy as R E(t, x) = 12 R3 mv2 f (t, x, v) dv. which is a macroscopic variable. It is not the temperature itself since in that case we need one more condition to have the relation: E = 3kT n where k is the Boltzmann constant. In order to generalize and summing up the results let us put a table in which the above passage are represented. Summing up we have three class of description linked to three type of equation:

1.3 Modeling Scales and Representation

13

• microscopic −→ ODE • kinetic −→ Lf = N f • macroscopic −→ PDE As we saw, a state variable for macroscopic scales is a local average on microscopic states and it is called macroscopic observable. Local averaging is a good approximation when a sufficiently small volume contains a number of elements large enough. Then, in order to use the macroscopic representation, we have to suppose this approximation, which is called continuous matter assumption, to be true. In this way we can define macroscopic state variables as limits on a volume that can be as small as we like; for example, density mass ∆m ∆x→0 ∆x

(1.1)

∆E[v] , ∆x→0 ∆x

(1.2)

ρ(t, x) = lim or local mean velocity

v(t, x) = lim

where E(v) is the average speed of particles.

14

1. Mathematical Models and Modeling Scales

Finally, using the mesoscopic scale allows us to consider a single particle, which is said test particle: this can not be identified with a real particle, but at the same time is a representant of each particle of the system.

1.4

Elementary Examples of Models at Different Scales

Microscopic Scale For this scale we could present any model of rigid-body dynamics in fact in these cases we consider the single objects and the model can be classify as discrete-dynamic so it will be stated with ODE and it is a typical microscopic model. We will describe a simple problem of rigid-body dynamics with a mass and an elastic wire. To simplify the model, we can think that the mass is concentrated in the center of mass P of the object and the wire attached directly to P. The mass m is constrained to translate along a horizontal line and is attached to an elastic wire. If we call x the distance of the mass from his resting position, the state could be ¶ µ dx u = u1 = x , u2 = dt From the Newton’s principles of classical mechanics and from a phenomenological model of the elastic behavior of the wire we obtain ¯ ¯(p−1) ¯ dx ¯ d2 x dx m 2 = F − kx − c ¯¯ ¯¯ dt dt dt where k is the elastic constant, F is a constant force directed along the horizontal axis and the last member represent the friction force that depends on p. We have non-linearity if p 6= 1 or if k = k(x). The values of parameters k, p, F and m are related to the specified problem. So we can write the following system of two first order equations that represent the evolution of u

 du1   = u2    dt      du2 = F − k u1 − c |u2 |p−1 u2 dt m m m

1.4 Elementary Examples of Models at Different Scales

The model is consistent, since there are two linearly independent equations corresponding to the two components of the state variable. To determine the solution of this kind of problem we have to know the initial state of the system ¡ ¢ u(t0 ) = u0 = u10 , u20 . that represents the initial position and velocity of the mass. At this point, we can rewrite this model in a dimensionless form. This is an important step for every models because it gives us a deeper understanding of the system because in this way the causes and the effects have the same order of magnitude and can be compared. If we introduce a critical length L and a critical time T that will be defined subsequently, we can make the following change of variable x ¯=

x ; L

t t¯ = . T

Substituting in the linear model we obtain d2 x ¯ kT 2 (1 − x ¯) = dt¯2 m and now it comes naturally to define L as the maximum length of the wire and T as the time that make the coefficient kT 2 /m = 1. p is L = F/k and T = m/k.

The right choice

The evolution equation of x ¯ becomes d2 x ¯ = (1 − x ¯) 2 ¯ dt that gives us the possibility to resolve many other similar problems with different values of the parameters. So the obtained equation is universal and it has a clear physical interpretation, i.e. it stands for every phenomenon of this kind. Macroscopic Scale As a simple model at macroscopic scale we will present a one-dimensional Heat Transfer in a rod of length L. The state variable u is the temperature that is an averaged quantity and that depends not only on time but also on

15

16

1. Mathematical Models and Modeling Scales

space so in this case we will have a dynamic-continuous model described by a PDE. u = u(t, x),

x ∈ [0, L]

If we introduce the heat flow q, we can write the balance equation cA

∂u ∂q dx = −A dx ∂t ∂x

where c is the heat capacity and A is the cross section of the rod. The phenomenological relation is q = −h

∂u ∂x

where h is the heat conduction coefficient. It comes that if c and h are constant, we have the following linear parabolic second-order PDE ∂u ∂2u = k 2, ∂t ∂x

k=

h = const c

For this type of equations we must know the boundary values u0 (x) = u(t0 , x), a(t) = u(t, 0),

∀x ∈ [0, L]

b(t) = u(t, L),

∀t > t0

Mesoscopic Scale Let us consider a community of people: we re interested in distribution of wealth among people. Let us define a dimensioneless microscopic variable u ∈ I = [0, 1].It corresponds to the amount of money owned by each subject: u = 0 means no money, while u = 1 represents the richest person. Let f (t, u) be the probability density obtained dividing the probability distribution of u by the number of subjects, then we have Z

u2

u1

f du = P(t, u ∈ [u1 , u2 ])

0 ≤ u1 ≤ u2 ≤ 1;

(1.3)

the mean value of u is Z E[u] =

uf du. I

(1.4)

1.5 Critical Analysis

17

Subjects, seen as groups of individuals, are continuously in competition: we can use a term η that measures the quickness or slowness of these competitions. In fact, these competitions take place through interactions among economic classes and η represents their frequency. Then let ϕ(u1 , u2 ; u) be the probability that, if the economic class u1 interacts with u2 , the subjects with economic condition u1 get to a new condition u. If we know η and ϕ, we can build a model starting from a balance equation within the elementary volume of the microscopic states du. Let G be the inlet flux of gain and L the outlet flux of loss: their difference leads to variations of f ∂f du = (G − L)du. ∂t

(1.5)

If we know ϕ, we can easily obtain G and L as Z G(t, u) = η I×I

ϕ(u1 , u2 ; u)f (t, u1 )f (t, u1 )du1 du2

(1.6)

Z

L(t, u) = f (t, u)η I

f (t, u1 )du2 ;

(1.7)

so, if we insert η in the dimensionless variable t, the final model is

∂f (t, u) = ∂t

1.5 1.5.1

Z I×I

ϕ(u1 , u2 ; u)f (t, u1 )f (t, u1 )du1 du2 − f (t, u).

(1.8)

Critical Analysis Boundary Conditions

Generally a mathematical model is an evolution equation which describe the evolution of some selected aspects of the real system and it can be defined a problem only if it is linked to the so-called initial and boundary condition. Indeed even the simplest differential model cannot predict the future if the behavior on the boundaries and at initial time are unknown. So the description of the system is obtained solving the mathematical problem generated by the application of the model with the values of the parameters found after physical experiments. In Section 1.2 we gave the definition of mathematical

18

1. Mathematical Models and Modeling Scales

problem in which we spoke about initial or boundary conditions now let us give some remarks about them. Boundary Conditions for ODE Firstly, let us consider microscopic models defined in the time interval [0, T ]: ordinary differential equations like

dui dt

= f (t, ui ) for each object can

be joined to initial conditions, like ui (t = 0) = ui0 or to limit conditions, such as ui (t = T ) = uiT . Generaly if we consider u,f ∈
du = f(t, u) dt

As we have seen in the first example of 1.4, a problem stated with an nth order ODE is well-formulated if we have n values of u or of his derivatives because we have to determine n constants. However not all the possible values make the problem well-posed. In fact the well-posedness is strictly related also with the forcing function f. For example if f is a continuous Lipschitz function we can give the initial values of u and all his derivatives until the (n − 1)th order calculated in t = 0 to have the well-posedness. Some problem could give us more values of u (for example the values in t = 0 and in t = T ) instead of his derivatives. In these problems the Lipschitz condition for f is sufficient to assure the existence but not the uniqueness of the solution. Obviously all the linear models satisfy the Lipschitz condition so they are all well-posed. Boundary Conditions for PDE We have already seen in the second example of 1.4 that boundary conditions for PDE are more complicated because we have a more complex domain for the state variable so we need also the values of the state variable at the

1.5 Critical Analysis

19

boundary of the space domain. A PDE could have Dirichlet Boundary Conditions, Neumann Boundary Conditions or a linear combination of them (Robin Boundary Conditions). The former is a condition that involves the state variable, the latter involves his space derivatives. It’s difficult to give general conditions for the well-formulation but generally we can give the right boundary information studying the single equation and his characteristics curves. For example a parabolic and a hyperbolic equation have different correct boundary conditions. To assure the well-posedness we can use instead some important theorems like the Lax-Milgram Lemma that demands not only hypothesis on the boundary conditions but also on the structure of the differential operator (continuity and coercivity). If both time and space derivatives appear in the mathematical model, then both initial and boundary conditions usually have to be assigned. The relative mathematical problem is then called initial-boundary-value problem. If, instead, the mathematical model is time independent no initial conditions are needed. In this case, the mathematical problem is only defined boundary-value problem. Summing up we can say that some different types of initial-boundary conditions can be given to a model which can be written as ∂u ∂2u ∂2u ∂u ∂2u = f (t, x, , ..., 2 , 2 , ..., , ...) ∂t ∂x1 ∂x1 ∂x2 ∂x1 ∂x2 where u = u(t, x) : [0, T ] × D → < is the dependent variable and the boundary of the domain of D ⊆ <3 of the space variables is denoted by ∂D. • Dirichlet problem: the initial-boundary value problem for the scalar equation is stated with initial condition u(0, x) = ϕ(x), ∀x ∈ D, and boundary conditions ∀t ∈ [0, T ] , ∀x ∈ ∂D : u = α(t; x ∈ ∂D), given as functions consistent, for t = 0 with the initial condition above. • Neumann problem: the initial-boundary value problem for the scalar equation is stated with initial condition u(0, x) = ϕ(x), ∀x ∈ D,

20

1. Mathematical Models and Modeling Scales

and boundary conditions ∀t ∈ [0, T ] , ∀x ∈ ∂D :

∂u = γ(t; x ∈ ∂D), ∂n

given as functions of time and where n denotes the normal to ∂D directed outside D – Mixed problem: the initial-boundary value problem for the scalar equation is stated with the same initial condition and to Dirichlet boundary conditions on one side of the domain and Neumann ones on the other side, such as u(t, ∂D1 ) = α(t)

and

∂u (t, ∂D2 ) = γ(t)∀t ∈ [0, T ] ∂x

or ∂u (t, ∂D1 ) = δ(t) ∂x

and

u(t, ∂D1 ) = α(t)∀t ∈ [0, T ]

But not every combination of Dirichlet and Neumann condition gives a well-posed problem, for example, if we consider a onadimension problem: u(t, a) = α(t) ∀t ∈ [0, T ] ∂u ∂x (t, a)

= γ(t) ∀t ∈ [0, T ],

they are given on the open boundary, constituted by the point x = a, so they do not have a physical meaning, and they takes a not stable solution. • Robin problem: the mixed initial-boundary value problem for the equation is stated with the same initial condition , to the boundary conditions defined as a linear combination of Dirichlet and Neumann boundary conditions , such as c1 (t)α(t) + c2 (t)δ(t) = ua (t)

and

c3 (t)β(t) + c4 (t)γ(t) = ua (t)

where c1,2,3,4 and ua,b are given functions of time.

1.5 Critical Analysis

21

• Non linear boundary-value problem: in this case the boundary conditions are defined as a nonlinear combination of Dirichlet and Neumann boundary conditions, say g1 (α(t), δ(t)) = ba (t)

and

g2 (β(t), γ(t)) = bb (t)∀t ∈ [0, T ]

where ba,b are given functions of time and g1,2 are suitable functions of their arguments. Obviously we need a number of boundary conditions equal to the highest order of the partial derivative with respect to the dependent variable to have a well-formulated problem. Some problems may refer to systems in an unbounded domains, in these cases boundary conditions have to be stated as above at the boundaries x → ±∞ . Asymptotic behavior can be stated letting the time variable to infinity. If the model is defined by a system of equations, then initial and boundary conditions have to be assigned for each equation according to the rules stated above (also the rules about unbounded domain).

1.5.2

Complexity

Referring now to complexity problems in modeling, we can say that every mathematical model tries to reduce the complexity of the real system through a simplified description by a finite number of variables. So on one side it is obvious that increasing number of variables (generally linked to the complexity of real world) makes the model closer to reality but on the other it could take computational and complexity problems. Simple examples shows that the computational time increases exponentially (or explodes) with respect to the number of variables. This problem influences also the choice of the observation scale as we will see in the next chapter.

Chapter 2

Models at Macroscopic Scale 2.1

Introduction

As we have introduced in 1.3, macroscopic models are characterized by a state variable that depends also on the space so they need a different approach from the microscopic ones. In this chapter we will analyze these models and in particular the diffusion ones, but all the problems of continuum physics are modeled at this scale. Macroscopic models are based on the continuum matter assumption: it states that given two points belonging to the system, however close each other, some matter is included in the region between this two points. If the distances between the single elements of the system is sufficiently small with respect to the characteristic dimensions of the body, the continuum hypothesis is an acceptable simplification, otherwise it is better to deduce a model at microscopic scale. • In 2.2 we give an example of finding out a diffusion macroscopic model, and after we propose a general method to obtain an evolution equation and we analyze some conservation and equilibrium equations. • Section 2.3 deals with the description and the properties of diffusion models and different vehicular traffic models are presented. • Section 2.4 speaks about the classification of the different PDE models and the different type of initial-boundary conditions used to formulate the mathematical problem. 23

24

2. Models at Macroscopic Scale

2.2

Deduction of Models

In the construction of a mathematical model we can define five different steps: • Step 1: Select the state variable u. • Step 2: Model the interactions between the inner system and outer environment, for example forces, stress fields, fluxes, sources, etc. • Step 3: Postulate the equilibrium, conservation, and/or balance laws which involve the state variable. • Step 4: Select suitable phenomenological relations(or constitutive laws), that could model the material behavior. • Step 5: Derive a model which can describe the evolution of u using the equations of Steps 3 and 4. As example of deduction of the model we consider a duct filled with a fluid (for example, but we’d consider other type of matheria) at rest and a pollutant diffusing in the duct in the direction of x-axis,so that we have a one-dimension problem. We assume the following conditions: we call c(x, t) the concentration of the pollutant and ρ0 its the density, which we can consider constant. As told in step 3 we have to choice some balance equilibrium or conservation equation. We consider a mass conservation equation. Considering the flux q = q(t, x) ,let q − the outlet and q + the inlet fluxes of fluid. They can be ralated by the equation q+ = q− +

∂q ∂x

dx,

Now we can write a balance equation, equating the net flux rate trough the section A of the duct to the increase of mass in a small volume, i.e. A dx. Then we obtain the following equation: ∂q − + ρ0 A ∂c ∂t dx = A(q − q ) = −A ∂x dx.

2.2 Deduction of Models

25

Now we have to join the flux with the concentration to obtain a phenomenological (material) model, which is not linked to the structure of the material but to the phenomenon observed. In this case we use a Fourier law: ∂c , q = −h0 ∂x

where h0 is the diffusion coefficient, that at present we suppose to be constant. This law doesn’t have universal validity. Now, substituting the above equation in the conservation one, we obtain : ∂c ∂t

=

h0 ∂ 2 c ρ0 ∂x2 ,

that is a linear parabolic partial differential equation. We can obtain the dimensionless equation by introducing the following quantities and relations: • u=

c cmax

where cmax is the maximam concentration of the fluid;

• t=

tvero Tc ,

where Tc = ρ0 L2 /h0 is the critical time and L is the length

of the duct; • x=

xvero L .

So by easy passages we obtined: ∂u ∂t

=

∂2u ∂x2

A deeper analyze can be done for example if we considered h0 not constant,which is more closely to the reality, in fact the higher is the concentration, the harder is diffusion. Moreover, if the concentration is zero, there is no diffusion (it is not necessary), such as if the concentration is maximum. If we consider a nonhomogeneous material,we can introduce a diffusion coefficient k depends on the space variable: k = h0 h(x), which, calling u = c,leads to ∂u ∂t

=

∂ ∂x

µ ¶ ∂u h(x) ∂x ,

26

2. Models at Macroscopic Scale

that is still a linear equation. If we want to consider a model that also deals with the situations of null concentration or saturated material, we we have to make k dependent on the concentration u: k = h0 h(u). An example can be: h(u) = uα (1 − u)β , with α, β ≥1. If α=β=1, we obtain a cubic nonlinearity in the equation: µ ¶2 ∂u ∂2u ∂u . ∂t = u(1 − u) ∂x2 + (1 − 2u) ∂x If our system is far from the equilibrium state or if there is a big variations of u, for example a termic shock in a heat transfert model we should consiteder that h may depend also on the gradient ∂u/∂x: µ ¶ ∂u k = h0 h(u) + γ . (2.1) ∂x But, when h → 0, the effect of the gradient is not so important; so we can correct the preceding expression: µ ¶ ∂u k = h0 h(u) 1 + γ . (2.2) ∂x Of course, if we choose this expression for k, we will have to assign one more boundary condition. It is also important to say that different causes produce different effects. In order to see this aspect let us consider the heat diffusion model in a bar: it has the same form seen until now, but here u represents the temperature. In our model the cause always produces an increase or a decrease of temperature. Another important consideration is that we may have an acceleration: ∂u ∂t

+

2 ε ∂∂t2u

=

∂ ∂x

µ ¶ ∂u h(u) ∂x ,

where ε is a parameter which measures the attitude of the material to change its temperature.

2.3 General deduction and different type of flux

2.3

27

General deduction and different type of flux

As we have seen in chapter 1 macroscopic models come from a stochastic analysis and an averaging of the microscopic spatial model. Generally we can define the position x and the velocity v of the elements of the system. If we accept the continuum assumption we can integrate the function f (x, v, t), that is the probability density function and we can obtain the following macroscopic observables Z ρ(x, t) = f (x, v, t)dv

Z q(x, t) =

vf (x, v, t)dv

where ρ and can be interpreted as a density (of mass, energy, number of individuals, etc.) and q is the flux. At this point we can postulate the mass balance equation. It can be stated like d dt

Z

Z ρ(x, t)dV =

P

Z q(x, t) · ndσ +

∂P

r(x, t)dV P

where P ⊂ Ω is an arbitrary volume and r is a ”source” term; if this term is null we have a mass conservation. From this, using the Reynolds Transport Theorem and the Gauss Divergence Theorem, we obtain ∂ρ +∇·q=r ∂t

(2.3)

This is the first important equation that is valid for most models. For ¡ ¢ many macroscopic models the state variable is u = ρ(x, t), v(x, t) where v = v(x, t) is the mean velocity, so the mass conservation can also be stated as

∂u + ∇ · ρv = r ∂t Sometimes (for example in continuum mechanics and fluid dynamics)

there is another important equations that comes from the assumption of the conservation of momentum. However, if this assumption is not valid, is sufficient to consider an operator F = F(x, v, t) that is a phenomenological information and that have to be determined for each particular system. So this equation can be stated as Z Z Z d qdV + ∇(q · v)dV = FdV dt Ω Ω Ω

28

2. Models at Macroscopic Scale

and with the theorems cited above we could write ∂ (q) + ∇(q · v) = F ∂t finally with the conservation of mass this equation simply reduces to ∂ ˜ (v) + v · (∇ · v) = F ∂t If we can’t postulate the conservation of momentum the model is not autosufficient and we have to write another equation. Other relations we can use in Step 3 could come from the Energy balance, from other experimental observations (for example the Coulomb law in Electromagnetism) or from modelizations of the microscopic interactions in the system which can give us information on the macroscopic model (for example in rarefied gasdynamics, porous media or population dynamics). All these relations with the constitutive laws can give us the possibility to derive many different models. Suppose we have to study the evolution of a density u. In the previous chapter we called it ρ. If we can’t apply the conservation of momentum, the balance equation 2.3 still not represent directly a differential equation on the variable u. To ”close” the problem is necessary to describe in other ways the flux q, by modeling the so-called phenomenological relations (or constitutive laws). If we consider also the source term r we must find relations of this type: q = q(u, ∇u)

r = rint (u, ∇u) + rext

where rint is the internal production of mass (or individuals), for example caused by chemical reactions or birth and death of populations. The term rext doesn’t depend on u and corresponds to an external production (for example immigration in the population models). Using these relations the mass balance becomes a PDE on the only u: ∂u + ∇ · q(u, ∇u) = rint (u, ∇u) + rext . ∂t

(2.4)

At this point, by neglecting the source term, the choice of the phenomenological relation on q is very important . We report some examples:

2.3 General deduction and different type of flux

• Linear flux q = uc, where c = c(x, t) is an assigned function called transport velocity that represent the velocity of the particles in each point. In this case we obtain a transport equation. • Convective flux q = g(u). In this case the velocity doesn’t depend linearly on u and we obtain a convection-transport equation • Diffusive linear flux q = −k∇u, where k is a constant called diffusivity. This relations describe the diffusion phenomenons. We can see that the flux is orthogonal to the level curves u = const and the vector q is directed where the function u decreases. We have a mass transportation from high density to low density places. The resulting equation will be ∂u = k∆u ∂t where ∆x is the Laplacian operator. • Diffusive non-linear flux q = −k(u)∇u. This time the diffusivity k depends on u. The diffusion equation can be written as ¡ ¢ ∂u = ∇ · k(u)∇u . ∂t Instead if k depends only on x we still have a linear equation but we will have one more term depending on (∇u)2 .

29

30

2. Models at Macroscopic Scale

2.3.1

Vehicular Traffic Model

Now let us give an example that shows how the same physical system can be modeled with different models. Let us consider the one dimensional flow of vehicles along a road of length L which is very long and in which travelling vehicles can change lane. We want to model it at a macroscopic scale so we must suppose that the distance between the vehicles is small and that the road is sufficiently long so that the continuum matter assumption is a good approximation. As we have seen in section 2.2 we can consider a density ρ and a mean velocity vt .The diffusion and transport model is ∂ρ ∂tt

+

∂ ∂xt (ρvt )

= 0.

Introducing new dimensionless analysis u=

ρ ρmax

;

u=

vt ; vmax

x=

xt ; L

t=

tt ; T

we can write the balance law as ∂u ∂ + (uv) = 0. ∂t ∂x Now we need a phenomenological relation, such as v = v[u]. We can initially consider a model in which v = 1 − u so we obtain this equation: ∂u ∂u + (1 − 2u) = 0. ∂t ∂x Comparing this model with the reality we can find two important defects: first we have that if u → 1/2, then ∂u/∂x → ∞. This means that there is an unacceptable discontinuity and this does not reflect the reality. The second problem is that the model formally needs only one boundary condition, but we would have a boundary condition at the beginning and one at the end of the road. To modify our model we can observe that, if we suppose that the road is sufficiently wide, there is a diffusion phenomenon, so we can add a non linear diffusion term and the equation becomes: µ ¶ ∂u ∂ ∂u ∂u = (2u − 1) +ε h(u) . ∂t ∂x ∂x ∂x Now the solution is smooth enough but the diffusion term that we added is not justifiable and we can’t give a method to calculate ε and h(u)

2.3 General deduction and different type of flux

Another important traffic model is the De Angelis one. It comes from the observation that a driver chooses its own velocity looking at the medium density that he can only see in the part of road in front of himself and not all around. So the driver is influenced by the apparent density µ ¶ ∂u u∗ = u 1 + ε(1 − u) . ∂x From this equation we can see that f u is next to 0, the model brings to a very high diffusion, independently on the value of the gradient. Using the phenomenological relation v = 1 − u∗ , the De Angelis model lets us write an equation of the type: µ ¶ ∂u ∂u ∂ ∂u = (2u − 1) + k(u) , ∂t ∂x ∂x ∂x where k(u) = εu2 (1−u) Now we can suggest a method to calculate ε,starting from the difference between v in stationary uniform conditions and v in stationary non-uniform conditions. Finally we can observe that De Angelis model still provides infinite propagation velocity. For this problem we can find these solutions. • Choose k(u) to obtain a degenerate parabolic equation with finite propagation velocity. • Add to the model an acceleration term so that it could be written the following hyperbolic equation µ ¶ ∂u ∂u ∂ ∂u = (2u − 1) + k(u) , ∂t ∂x ∂x ∂x . But the coefficient µ is not experimentally calculable. • Use, with the mass conservation equation, also the conservation of momentum to obtain  ∂u ∂   + (uv) = 0   ∂x  ∂t      ∂v + v ∂v = F [u, v] ∂t ∂x

31

32

2. Models at Macroscopic Scale

2.4

Classification and Formulation and diffusion model properties

Once a mathematical model has been deduced , the statement of the mathematical problems (generated by the application of the model to the real world) requires dealing with a detailed analysis of the qualitative properties of the model. Indeed also the behavior of the solutions and the application of computational schemes depend on the properties of the model which can be classified according to its mathematical structure. As we have seen in this chapter a large variety of macroscopic models can be written in terms of partial differential equations. Consider the quite general linear second order scalar partial differential equation in n independent variables, written as follows : n X i,=1

n

X ∂u ∂2u (x) + Ai (x) (x) + a0 (x)u(x) = F (x) Aij (x) ∂xi ∂xj ∂xi i=1

which can be formally written as : Lu(x) = F (x) where the dependent variable x may include time. To simplify the notation let us introduce the square matrix A = (Aij )1≤i,j≤n which is symmetric , since

∂2u ∂2u = ∂xi ∂xj ∂xj ∂xi

(and then Aij = Aji )) so that her eigenvalues are real. Then above equation can be rewritten as: ∇x · (A∇x u(x)) + ∇x · (au(x))) + a0 u = F (x) with a = (a1 , a2 , ..., am )T , the coefficient of the first order part. The classification is proposed according to the structure of the principal part , the terms of the equation involving the higher order derivatives i.e., of the partial differential equation. This classification is obtained analyzing the eigenvalues of the matrix A. In particular we analyze the following three important situation.

2.4 Classification and Formulation and diffusion model properties 33

• All the eigenvalues of A are different from zero and one eigenvalue has a different sign respect the others; in this case the operator L is called hyperbolic. • One eigenvalue of A is zero while the others have a constant sign; in this case the operator L is said parabolic. • All the eigenvalues of A are different from zero and they have the same sign; the operator L is called elliptic. Summing up we have: ∀i, λi 6= 0, ∃!λi > 0 ∀i, λi 6= 0, ∃!λi < 0 ∃λi = 0 ∀i, λi > 0 ∀i, λi < 0

) ⇒ hyperbolic equation, o ⇒ parabolic equation,

(2.5)

) ⇒ elliptic equation.

This classification influences the admissible kind of boundary and-or initial conditions for the equation, some properties of the solution and the techniques used to solve the problem. Now let us give some remarks about this classification: • if the coefficients of the matrix A depend on the space variable, the type of equation may change from point to point; • for non linear equations the type of degeneration may depend not only on the point in space but on the solution itself. However by a linearization of the equations we can extend the above classification. So we can say that a non linear operator is of the same type of the related linear one; • the above classification is exhaustive only in two dependent variables since the matrix A is 2 × 2 and it has almost two eigenvalues; we have to say that in this case it’ s possible another way to classify the equation, in fact we have: 2 A(t, x) ∂∂t2u

+

∂ 2u 2B(t, x) ∂t∂x

+

2 C(t, x) ∂∂xu2

µ ¶ ∂u ∂u = F t, x; u, ∂t , ∂x .

34

2. Models at Macroscopic Scale

and the classification is determined by the sign of the discriminant B 2 − 4AC : B 2 − 4AC > 0 ⇒ hyperbolic equation B 2 − 4AC = 0 ⇒ parabolic equation

(2.6)

B 2 − 4AC < 0 ⇒ elliptic equation. We can observe that the discriminant depends on independent variables, so its sign can change and then also the equations can change type. • the level curves in the (ξ1 , ξ2 ) plane of the associated quadratic form Q(ξ) = ξ T A ξ

with

ξ = (ξ1 ξ2 )T

are ellipses or degenerated parables or hyperboles depending by the type of the operator (since the name of the above classification). It can be shown that the qualitative behavior of the solutions depend on the properties of the model according to the above classification. Parabolic equations describe diffusion-like phenomenons. In general, the solution of a parabolic problem is smooth with respect to both space and time even if the initial data are not continuous. In other words parabolic system have a smoothing action and singularities do not develop nor can be maintained. So , even if the initial data have compact support, at any t > 0 the initial condition propagates everywhere. This is usually indicated by saying that parabolic systems have an infinite speed of propagation. Elliptic equations describe systems in equilibrium or steady state. The equations describe, for instance, the final state reached by a physical system described by a parabolic equation after the transient term has died out. Hyperbolic equations describe wave-like phenomena. The solution of hyperbolic initial-value problem cannot be smoother than the initial data.On the other hand , it can be develop, as time goes by, singularities even from smooth initial data, which then characterize the whole evolution and are propagated along special curves called characteristics. In particular if the solution has initially a compact support then such a support expands with a finite speed of propagation.

Chapter 3

Discretization Methods of Macroscopic Models 3.1

Introduction

This chapter deals with the theory of the discretization methods for continuous model. • In Section 3.2 we give an idea of the finite differences method, collocation method and domain decomposition methods. In particular we apply these methods to a parabolic equation (advection-diffusion). • Section 3.3 is formed by some examples of simulation of diffusion model. In order to do this we have used the software Mathematica which allows to simulate the models and to plot the graphics of the solutions.

3.2

Discretization of Diffusion Models

Once the mathematical problem has been correctly formulated and studied at the qualitative level, we can turn our attention to the numerical simulation. We have to say that there is a large variety of numerical methods dealing with partial differential equation. In this section we give an idea of three methods: the finite difference , the collocation and the domain decomposition one. 35

36

3. Discretization Methods of Macroscopic Models

3.2.1

The finite difference methods.

The first step of the finite-difference methods is to divide the interval [a, b] into N subintervals , usually of the same length ∆x, a = x0 ≤ x1 ≤ ... ≤ xN −1 ≤ xN = b The basic idea is then to replace every partial derivative that appears in the differential equation and in the initial and boundary condition by an appropriate approximation. Of course, this can be done in several ways. For instance, first order space derivatives can be approximated as u(t, x + ∆x) − u(t, x) ∂u (t, x) ∼ = ∂x ∆x or as

u(t, x) − u(t, x − ∆x) ∂u (t, x) ∼ = ∂x ∆x

Both of them are first order approximation called, respectively, forward and backward difference. Better accuracy is given by the central difference approximation u(t, x + ∆x) − u(t, x − ∆x) ∂u (t, x) ∼ , = ∂x 2∆x which is second order. As far as higher derivatives are concerned we recall the central differences approximations ∂2u 1 (t, x) ∼ [u(t, x + ∆x) − 2u(t, x) + u(t, x − ∆x)] = 2 ∂x ∆x2 ∂3u 1 (t, x) ∼ [u(t, x + 2∆x) − 2u(t, x + ∆x) + 2u(t, x − ∆x) − u(t, x + 2∆x)] = 3 ∂x ∆x3 ∂4u 1 (t, x) ∼ [u(t, x + 2∆x) − 4u(t, x + ∆x) + 6u(t, x) − 4u(t, x − ∆x) + u(t, x + 2∆x)] = 4 ∂x ∆x4 So the system of partial differential equation can be written as: du (t, xj ) = f (t, xj ; u(t, xj ), δx u(t, xj ), δxx u(t, xj ), ...) dt which is an ordinary differential equation approximating ∂u ∂u ∂ 2 u (t, x) = f (t, x; u, , , ...) ∂t ∂x ∂x2

3.2 Discretization of Diffusion Models

and where the finite differences δx u, δxx u, ... are chosen to approximate ∂2u ∂x2

37 ∂u ∂x ,

and so on. We have to remark that the above ordinary differential

equations may be solved using the classical numerical methods for discrete initial-value problems, such as: • Euler method, • Crank-Nicolson method, • Runge-Kutta methods, • Adam-Bashforth or Adam-Moulton methods , and so on. Let us give you a definition: if in f all quantities are evaluated at time ti then we have an explicit method , if the right-hand side of the equation above is evaluated at time ti+1 an implicit method and if some of the terms in the right-hand side are evaluated at time ti and some at time ti+1 an semi-implicit method. To approximate the partial derivative in time we can use for example the forward difference approximation u(t + ∆t, xj ) − u(t, xj ) ∂u (t, xj ) ∼ = ∂t ∆t and if f is computed explicitly, we have u(t + ∆t, xj ) = u(t, xj ) + ∆tf (t, xj ; u(t, xj ), δx u(t, xj ), δxx u(t, xj ), ...)) If instead the time derivative is approximated by the central difference approximation u(t + ∆t, xj ) − u(t − ∆t, xj ) ∂u (t, xj ) ∼ = ∂t 2∆t we have u(t + ∆t, xj ) = u(t − ∆t, xj ) + 2∆tf (t, xj ; u(t, xj ), δx u(t, xj ), δxx u(t, xj ), ...) A finite difference approximation is said to be nth order in space and mth order in time if approximations of at least order n are used to approximate space derivatives and approximations of at least order m are used to the time ones. The kind of approximation used strongly depends on the type

38

3. Discretization Methods of Macroscopic Models

of the equation considered. Indeed, the parabolic problems smooth out the solution, therefore the introduction of an artificial viscosity or the dispersion effect is not so important.The main trouble is that the stability requirements for explicit schemes are severe so that implicit schemes are preferred.Elliptic problems describe systems in equilibrium or steady states. Once every derivative is replaced by its finite-difference approximation, the partial differential equation is replaced by a system of equations in the unknowns u(xi , yj ) to be solved at the same time. Therefore the problem reduces to the solution of a large system of equations, This is usually done by iterative methods, such as Gauss-Seidel, Successive Over Relaxation (SOR) and so on. Instead hyperbolic problems govern propagation-like problems with possible formation of shocks and conservation of some physical quantities. So we have possibly to use a method that does’t smooth out discontinuities of the solutions or of its derivative, is able to preserve the physical quantities as required by the conservation equation and to detect properly the propagation speed. Another difficulty that occurs using finite difference methods is that in two or three-dimensional problems the space grid will not fit the domain, this implies a difficulty in properly giving the boundary conditions without loosing accuracy. To overcome this we may use suitable interpolation formulas to set boundary conditions or map the domain into a rectangle or use domain composition methods. Now let us give an example of using finite differences for parabolic equations. Consider the advection-diffusion equation ∂u ∂2u ∂u = ν(t, c) 2 − c(t, x) + f (t, x)u + g(t, x) ∂t ∂x ∂x

x ∈ [a, b]

with initial condition such as u(0, x) = uin (x) and Dirichlet boundary conditions: u(t, a) = ua (t)

and

u(t, b) = ub (t).

Using central differences for the space derivatives we can write: uj+1 − 2uj + uj−1 uj+1 − uj−1 duj = νj (t) − cj (t) + fj (t)uj + gj (t) dt ∆x2 2∆x

3.2 Discretization of Diffusion Models

39

for j = 1, ..., N −1 ; where now uj (t) = u(t, xj ). The equation above can be written as µ ¶ µ ¶ µ ¶ duj cj νj 2νj cj νj = + uj−1 + fj − uj + − + uj+1 . dt 2∆x 2∆x2 ∆x2 2∆x 2∆x2 So we have approximated a partial differential equation with a system of ordinary differential equation (which can also be written in vector form). Now we have only to choose a stable and sufficient accurate ODE solver.

3.2.2

The collocation method.

As for the finite difference methods, the first step of a collocation method is to select N +1 nodes in the space interval [a, b] a = x0 ≤ x1 ≤ ... ≤ xN −1 ≤ xN = b. Then the second step is to select a set of cardinal basis functions φj (xi ) such as the Lagrangian or the trigonometric polynomials or the fundamental splines so that φj (xi ) = δij

i, j = 0, ..., N

where δij is the Kronecker delta. The solution u(t, x) to the initial-boundaryvalue problem can then be approximated as def u(t, x) ∼ = uN (t, x) =

N X

uj (t)φj (x)

j=0

where uj (t) = u(t, xj ). The partial derivatives of u can be naturally approximated by the partial derivative of uN N

X duj ∂uN ∂u (t, x) ∼ (t, x) = (t)φj (x) = ∂t ∂t dt j=0 N

X dφj ∂u ∂uN (t, x) ∼ (t, x) = uj (t) (x) = ∂x ∂x dx j=0 N

X d2 φj ∂2u ∂ 2 uN ∼ u (t) (t, x) (t, x) = (x) = j ∂x2 ∂x2 dx2 j=0

40

3. Discretization Methods of Macroscopic Models

Substituting the above formulas into the sample advection-diffusion model ∂u ∂u ∂2u + c(t, x; u) = ν(t, x; u) 2 + F (t, x; u) ∂t ∂x ∂x we have: N N N N N X X X ¡ ¢ ¡ ¢X ¡ ¢X dφj d2 φj duj uj uj +F t, x; u φ φj +c t, x; uj φj = ν t, x; uj φj j j dt dx dx2 j=0

j=0

j=0

j=0

j=0

The second step is to consider the above equation in the collocations points xi so they can be written as: N

N

j=0

j=0

X X dui (t)+c(t, xi ; ui (t)) Dij uj (t) = ν(t, xi ; ui (t)) Bij uj (t)+F (t, xi ; ui (t)) dt where Dij =

dφj (xi ) dx

Bij =

d2 φj (xi ) dx2

are the components of two matrices D and B that can be computed once and for all as the collocation points xi and the cardinal basis functions φj are chosen. There are now two choices left arbitrary : • the collocation points; • the cardinal basis functions. Naturally the choice depends on the properties of the solution and on the accuracy with which the derivatives are computed. From former numerical computation courses we know that a Chebyshev distribution points is the best choice to couple with Lagrange polynomials, instead if we know that the solution is periodic, the natural choice is to use trigonometric polynomials with equidistant nodes. The former case deals with the Chebyshev collocation method, the latter as the Fourier collocation method. Before going on into the analysis oh these methods let us give another definition. Spectral accuracy: the adjective spectral is used to indicate an approximation with an error that decreases faster than any power of the number of collocation points for infinitely differentiable functions.

3.2 Discretization of Diffusion Models

41

Fourier Collocation Method In Fourier collocation methods the reference interval is [0, 2π] , so we have to map our interval [a, b] into it with a smooth transformation ξ = f (x) as: ξ = 2π x−a b−a or x = a +

(b−a) 2π ξ

It is possible to prove that the choice of

equidistant collocation points ξh =

hπ N

h = 0, ..., 2N − 1

is equivalent to approximating u as N −1 X

u(t, ξ) ∼ = uN (t, ξ) =

u ˜k (t)eikξ

k=−N

where i is the imaginary unit and u ˜k (t) =

2N −1 −ihkπ 1 X u(t, ξh )e N 2N

k = −N, ..., N − 1

h=0

So we obtain : N −1 X

u(t, ξh ) ∼ = uN (t, ξh ) =

u ˜k (t)e

ihkπ N

k=−N

for h = 0, ..., 2N − 1. These two transformation are usually called Discrete Fourier Transform (DFT) and Inverse Discrete Fourier Transform (IDFT). The values of the derivatives in the collocation points can be computed from the values of u in the collocation points by the following steps: compute u ˜k knowing u(t, ξh ) and than derive to obtain N −1 X ∂uN ik˜ uk (t)eikξ (t, ξ) = ∂ξ k=−N

(1)

which says that the coefficients u ˜k of the derivative are: (1) u ˜k

2N −1 −ihkπ ik X = ik˜ uk = u(t, ξh )e N 2N h=0

for k = −N, ..., N − 1. At the end we have to transform back to finally obtain:

N −1 X ∂uN (1) ihkπ (t, ξh ) = u ˜k e N ∂ξ k=−N

h = 0, ..., 2N − 1

42

3. Discretization Methods of Macroscopic Models

This approximation of the derivative is called Fourier collocation derivative and is spectrally accurate. The formula can be easily generalized to higher order derivative by: (n)

u ˜k = (ik)n u ˜k =

2N −1 −ihkπ (ik)n X u(t, ξh )e N 2N h=0

The same result can be achieved by matrix multiplication. In fact we can write:

2N −1 X ∂uN (t, ξh ) = Dhk uk (t), ∂ξ k=0

where Dhk =

(−1)h+k

if h 6= k;

2 sin (h−k)π 2N

or in vector form

else Dhk = 0

∂u (t) ∼ = Du(t) ∂ξ

where u(t) = {u0 (t), ..., uN (t)} is the unknown vector of the state variable n o in the collocation points and ∂u (t) ∼ = ∂u (t, ξ0 , ..., ∂u (t, ξN ) is the value of ∂ξ

∂ξ

∂ξ

the derivative of u in the collocation points. However despite its simplicity, this procedure is computationally longer. For high order derivatives ∂nu ∂ n uN ∼ (t) (t) = Dn u = ∂ξ n ∂ξ n where Dn is the nth power of the matrix D. Chebyshev Collocation Method In Chebyshev collocation methods, the reference interval is [−1, 1]. As before, the first step is to map the interval [a, b] into it with a smooth transformation, say the linear one ξ =1−2 or

x−a b−a

b−a (ξ − 1) 2 We remark that Chebyshev polynomials are characterized by the followx=a−

ing features:

3.2 Discretization of Diffusion Models

43

• range: [−1, 1]; • polynomials: T0 = 1; • weight: i = j 6= 0 :

T1 = x;

2/π ; 1−x2

Tm = 2xTm−1 − Tm−2 ;

i=j=0:

1/π . 1−x2

Instead Lagrange polynomials are given by the following expression: Li (x) =

(x − x0 )...(x − xi−1 )(x − xi+1 )...(x − xn ) (xi − x0 )...(xi − xn )

It is possible to prove that the choice of Chebyshev collocation points ξh = cos

hπ N

h = 0, ..., N

joined with Lagrange polynomials is equivalent to approximating u as def u(t, ξ) ∼ = uN (t, ξ) =

N X

u ˜k (t)Tk (x)

k=0

where Tk (x) are the Chebyshev polynomials and the coefficients u ˜k are given by N hkπ 2 X 1 u(t, ξh ) cos u ˜k (t) = N c¯k c¯k N

k = 0, ..., N

k=0

where c¯h = 2 if h = 0 c¯h = 1

if 1 ≤ h ≤ N − 1.

Not only this approximation is spectrally accurate, but is also computationally convenient since N

u(t, ξh ) = u (t, ξh ) =

N X k=0

u ˜k (t) cos

hkπ N

h = 0, ..., N

which enables the use of FFT algorithms in computing both the Discrete Chebyshev Transform (DCT) and the Inverse Discrete Chebyshev Transform (IDCT) and therefore in the evaluation of the derivatives. As for the Fourier collocation derivatives, the Chebyshev ones can be evaluated from the values of u in a collocation points using a scheme like this: knowing

44

3. Discretization Methods of Macroscopic Models

u(t, ξh ) to compute u ˜k ; derive the approximation to obtain the recurrence formula  (1)   u ˜N +1 = 0      (1) u ˜N = 0 (1) (1)   u ˜k = u ˜k+2 + 2(k + 1)˜ uk+1      (1) (1) u ˜0 = 12 u ˜2 + u ˜N +1

for k = N − 1, ..., 1,

Finally, transform back with the inverse Chebyshev transform, to obtain N

X (1) ∂uN hkπ (t, ξh ) = u ˜k cos ∂ξ N

h = 0, ..., N

k=0

This approximation of the derivative is called Chebyshev collocation derivative and is spectrally accurate. In general, for the nth derivative one has to repeat this scheme n times. It is useful to observe that the Chebyshev collocation derivative can also be computed using the matrix D defined above. As for the finite differences method we analyze the application of Chebyshev collocation method to the advection-diffusion model.  ∂u ∂u ∂2u   + c(t, x) − ν = F (t, x)   ∂t ∂x ∂x2    u(0, x) = u (x) in

  α ∂u   ∂x (t, a) + βu(t, a) = ua (t)    γ ∂u (t, b) + δu(t, b) = u (t) b ∂x We supposed ν constant. The collocation method leads to the semi-discrete equation  N X  duh  ˜ 2 uk = Fh  + ch D  hk   dt k=0 PN ˜  α  k=0 D0k uk + βu0 = ua    P  ˜ γ N k=0 DN k uk + δuN = ub

for

h = 1, ..., N − 1

where uk are the values of u in the collocation points µ ¶ kπ b−a cos −1 xk = a − 2 N

(3.1)

3.2 Discretization of Diffusion Models ˜ 2 is the second derivative matrix, which is the square of the first and D hk ˜ hk . Of course the initial-value problem has to be solved derivative matrix D with the initial condition uh (0) = uin (xh ). Obviously before choosing the ODE solver, we have to examine the eigenvalues related to the differential operators and to set the time step ∆t so that all them are squeezed into the stability region. Domain Decomposition Method The modeling and mathematical methods presented in the preceding sections were developed by adopting the same model and the same solution techniques over the whole domain of the dependent variable. However , this situation may be unrealistic in several physical situations. In fact, we have already seen that the partial differential equation(and even the mathematical model itself) may change type in different part of the domain. So the domain D might have a shape that does not possess the sufficiently regularity required for the application of the numerical scheme. Then a suitable subdivision of D may be organized in order to obtain a convenient applicability of the mathematical method that is adopted. For example collocation methods work well in rectangular domains or in domains that can be easily mapped into rectangles. So in general in this method we cut the domain D into m smaller pieces Di with: D=

m [

Di

i=1

averaging the state variable within each volume. It follows that, in each finite volume the dependent variable is ui = ui (t), and that the mathematical model is a finite one and consists in a set of n equations providing the evolution of the ui . So in each part of the domain we have a subproblem formed by the same conservation and-or equilibrium equations of the continuous model. The above procedure generates a finite model and the mathematical problem is stated in terms of a system of ordinary differential equations.

45

46

3. Discretization Methods of Macroscopic Models

3.3

Simulations of Diffusion Models

We have elaborated the following code in Mathematica that gives us the possibility of simulate a linear Diffusion problem with mixed boundary condition by using Chebychev collocation and Lagrange Polynomials as described in the previous section. In particular we used the equations 3.1. For each specified problem we have to insert the domain, the coefficients, the boundary and initial condition.

Mono Dimensional Heat Transfer The first simulation is taken from the Section 1.4. We want to simulate the heat diffusion in a rod of diffusivity k = 1, initially cold, heated quickly from one side. We expect that the bar starts heating from this side towards the other. We have a not flux condition on the not heated initially end (null Neumann initially condition) and a Dirichlet condition in the other

3.3 Simulations of Diffusion Models

side which is represented by a function which quickly increases and than becomes constant. Substituting the initial and boundary conditions, we obtain the following problem  ∂2u ∂u   =   ∂x2   ∂t u(x, 0) = 10   ux (0, t) = 0    √  17 u(1, t) = 30 t, with t ∈ [0, 1] and x ∈ [0, 1]. From the graphic of the solution we can see the result of the diffusion phenomenon.

The diffusion operator so tends to make the temperature of the not heated end equal to the temperature of the other end. The time needed to reach an equilibrium state depends on the diffusion coefficient. If we increase it, the time’d decrease. De Angelis model We have also simulated two problems of vehicular traffic by using the De Angelis model described in Section 2.3. We consider ε = 0.3.We know that

47

48

3. Discretization Methods of Macroscopic Models

trere is both transport and diffusion. The first simulation have an initial condition that have a maximum of density in the middle of the road and have constant boundary conditions, which means that we are modelling the rush-hour traffic. Obviously we should assume that it takes the same time to each toll collector serving a driver. The problem can be stated as  ∂u 1 2 ∂2u 1 ∂ 2 ∂u ∂u   = (2u − 1) + (u − u3 ) 2 + (u − u3 )   ∂t ∂x 3 ∂x 3 ∂x ∂x   2 −40(x−1/2)2 u(x, 0) = 0.2 + 5 e   u(0, t) = 0.2     u(1, t) = 0.2 with t ∈ [0, 1] and x ∈ [0, 1], which are the dimensioneless variables. From the plot of the solution we can see that the density (u(x, t) ∈ [0, 1] )tends to an uniform value in the whole road and the maximum moves and decreases quickly. The peak of the density moves in the direction of the vehicles. We can also observe that the big concentration of cars moves towards the going out toolbooth.

The second simulation of this model reflects the following situation: the

3.3 Simulations of Diffusion Models

initial density is the same of the fprmer simulation but while one boundary condition is constant, the other one has a sinusoidal behavior.  ∂u 1 2 ∂2u 1 ∂ 2 ∂u ∂u   = (2u − 1) + (u − u3 ) 2 + (u − u3 )   ∂t ∂x 3 ∂x 3 ∂x ∂x   u(x, 0) = 0.1   u(0, t) = 0.1     1 sin(10t) u(1, t) = 0.1 + 10 This plot gives us a good approximation of the reality in fact the sinusoidal

behavior is transmitted to all the road.At the end, we can observe alternatively stretches of road with high concentration of vehicles and stretches with a small density of cars, according to a sinusoidal trend of density depending on the spatial position. Intuitively we can say that this reflects what actually would happen.

49

Chapter 4

Conclusions Mathematical modeling constantly supports the development of applied sciences (for example biology, medicine end earth sciences) with the essential contribution of mathematical methods. Modeling processes are composed by well defined methods so that it is correct and useful to talk about the science of mathematical modeling. These complex procedures can be summarized in the following flow chart.

Analyzing it we can say that the first stage of the modeling process is the observation of the physical system which has to be modeled. With the term observation we also means organization of experiments suitable to pro51

52

4. Conclusions

vide quantitative information (values of parameters, for example) on the real system. The modelization process has to be considered as a sequence of operative stages which can be interrupted only when a satisfactory agreement is reached between simulations and observation of the phenomena. However we can point out some characteristic features of the modeling process: • generally a mathematical model, although approximating the physical reality should not hide relevant features , and in particular nonlinear behaviors. • Physical systems sometimes have stochastic behaviors. In some situations even if the mathematical model is of a deterministic type, the related mathematical problem may be of a stochastic type. For example if initial or boundary conditions or some parameters cannot be measured precisely we can give them affected by some stochastic noise. In those cases suitable mathematical methods need to be developed in order to deal with stochastic problems. • Modeling not only leads to a simulation of physical reality, but can also contribute to a deeper understanding of physical systems so that after the simulation the experimental observation can be revisited and hopefully improved. This work clearly shows that the science of mathematical modeling is deeply linked to the study of mathematical methods of functions approximation. So we can say that mathematical modeling requires observation, initiative and invention on one hand and knowledge theory of mathematical methods on the other. We try to put the attention also on the simulation which is an important tool in fact let us have the possibility to study how the reality evolves and changes under different hypotheses. Our proposal in this dissertation is not to give an exhaustive analysis of mathematical modeling science but to give some lines of conduct to follow . We have done this also throughout some examples. It is useful to stress that this is not a course of mathematical methods so sometimes we cite a method without studying it in detail, in these case we refer to [3] or other bibliography. So, for another time, it is emphasized the

53

close and deep connection between the mathematics (with his theory) and the real world which have to be analyzed and study by the applied sciences using tools such as mathematical model.

Bibliography [1] Bellomo N. and Preziosi L., Modelling Mathematical Methods and Scientific Computation, CRC Press, 1995. [2] Bellomo N., Da Angelis E. and Delitala M., Equations of Mathematical Physics and Applied Sciences, Birkhauser, 2005. [3] Quarteroni A., Modellistica Numerica per Problemi Differenziali, Springer, 2003.

55