Differential Equations Classification and Simulation November 25, - - PDF document

Differential Equations Classification and Simulation November 25, 2008

• M. Emmerich, LIACS

1

Classification

• Ordinary Differential Equations
• Systems of Coupled Ordinary Differential

Equations

• Partial Differential Equations

2

Ordinary differential equations

• Consider the equation:

dh dt = Q(t)/A(h) (1)

• Here h(t) describes the height of a fluid in a

reservoir depending on time, Q(t) describes the inflow in m3/s, and A(h) describes the area of the reservoir, which depends on the height.

• The equation describes an ordinary differ-

ential equation (ODE)

• In this lecture we will review methods to

solve ODEs numerically

3

Ordinary differential equation (2)

• An ordinary differential equation is an equa-

tion that involves one or more derivatives

• f an unknown function
• A solution of an differential equation is a

specific function that satisfies the equa- tion, e.g.

Equation Solution x′ = x x(t) = c et x′ = a x(t) = at + c x′′ + 9x = 0 x(t)c1 sin(3t) + c2 cos(3t)

• the letter c dentes an arbitrary constant
• this indicates that, in general, ODEs have

not a unique solution function

4

Initial-value problem

• Given an initial value x(t0) at some point t0

the solution curve of the ODE gets unique, e.g. Equation: x′ = x Solution: x(t) = cet Initial-value: x(2) = 5 Solution: x(t) = 5 e2et

• ODE problems with an initial value given

will be called initial-value problems

• Given the initial value, all past and future

states of the system described as ODE are determined

5

Boundary value problem −(y′)2 − 2by + 2yy′′ = 0 y(0) = a0, y(1) = a1

• Due to the type of condition, the problem

is a so-called boundary value problem

• The differential equation is of 2nd order,

as it includes first and second derivatives

• So-called shooting methods are used for

the numerical solutions of boundary value problems for ODEs

6

Solution of ODE

• Analytically, e.g.

dV/dt = at + b ⇔ (2) dV = (at + b)dt ⇔ (3)

• dV

=

• (at + b)dt

(4) V (t) = a 2t2 + bt + C (5)

• Symbolic solvers (Maple, DeRive, Mathe-

matica)

• de := D(x)(t) = 2 t + 4;

dsolve({ode, x(2)=5},x(t));

• Numerically (Euler, Runge-Kutta)

7

Coupled Differential Equations

• Example are Lotka Volterra Equations de-

scribing a predator prey system: x′ = ayx − b (6) y′ = −cyx + d (7) x is the population size of predator y is the population size of prey a: birth rate of predator (depends on y) b: death rate of predator c: death rate of prey (depends on x) d: birth rate of prey

• The system results in two solution curves

8

Other examples of Coupled ODE

• chemical reaction systems

state variables: concentrations of compo- nents constants:kinetic constants

• environmental models

state variables: CO2, temperature constant: energy intake from sun, heat absorption rates

• market models

state variables: price and demand constant: available ressources

10

Partial derivatives

• Consider a landscape described by the graph
• f function f(x1, x2)
• ∂f

∂x1 denotes the slope of the landscape in

direction

1

• ∂f

∂x2 denotes the slope of the landscape in

direction

1

• for d-dimensional function f(x1, . . . , xd), ∂f(x)

∂xi

denotes the slope of the landscape in the direction of the i−th unit vector ei

• for differentiable functions f:

∂f(x) ∂xi ( x0) = lim

h→0

f( x0) + f( x0 + h ei) h (8)

11

Gradient

• grad f(

x) = ( ∂f

∂x1, . . . , ∂f ∂xd) is the gradient

• f f at

x

• The gradient points in the direction of steep-

est ascent of f( x)

• The length of the gradient depends on the

steepness in that direction

12

Partial differential equations

• A partial differential equation is an equa-

tion with partial derivatives and typically also a derivative in time

• A very simple example:

(∂f ∂x, ∂f ∂y, ∂f ∂t )T = gradf = (sx, sy, st)T (9)

• Solution of this equation is a moving plane:

f(x, y, t) = sx · x + syy + stt + c (10)

• It has a well defined solution, if we provide

a condition f(x0, y0, t0) = h0

13

Partial differential equations

• Partial differential equations can also be

defined on vector fields

• Navier Stokes Equation for fluid flow of

compressible fluids are important example

• A closed form solution of these equations

is one of the millenium problems in math- ematics

• However, we can look at special cases, e.g.

incompressible fluids

• Or, we can solve the equations numeri-

cally using finite difference or finite ele- ments methods ∗

∗in turbulent regime difficulties will arise, however

14

Wave equation with u being the displacement

• f a particle:

∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 − ∂2u ∂t2 = 0 Heat equation with u being the local temper- ature: ∂2u ∂x2 + ∂2u ∂y2 + ∂u ∂z2 − ∂u ∂z = 0 Laplace’s equation, where u represents the ve- locity potential of a incompressible and irrotat- able fluid: ∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 = 0

15

Numerical solution of ODE

• Taylor methods can approximate a solution

to a differentiable function

• Differentiable functions can be locally ap-

proximated by the Taylor series x(t+h) = x(t)+hx′(t)+h2 2 x′′(t)+. . .+hn n!hnx(n)(t)

• when only terms up to order n are consid-

ered the method is called a Taylor method

• f order n

16

Euler’s method

• The Taylor method of order 1 is called the

Euler method

• To find an approximate solution to the ini-

tial value problem x′ = f(t, x(t)), x(a) = xa (11)

• ver the interval [a, b], the first two terms
• f the Taylor series are used:

x(t + h) = x(t) + hf(t, x(t)) (12)

• Formally this process is an iterated func-

tion system, and we could compute its Ljapunov coefficient

17

Pseudo code of Euler method

• The following pseudocode produces an ap-

proximation of the solution trajectory over the interval [a, b] with n points Euler method input: n, a, b, xa h ← (b − a)/n t ← a print 0, t, x for all k ∈ (1, . . . , n) do x ← x + hf(t, x) t ← t + h print k, t, x end for

• Eulers method is not very accurate.

The truncation error is O(h2) in one iteration, i.e. if we decrease h the error will decrease at a rate proportional to h2

18

Taylor method of order n

• Given x′ we first have to determine deriva-

tive functions x′′, . . ., x(n)

• Then we can use the following algorithm

to determine approximation over interval [a, b] Taylor method of order n input: a, b, xa, n h ← (b − a)/n t ← a print 0,t,x for all k ∈ (1, . . . , n) do x ← x + hx′ + h2

2 x′′ + . . . + hn n!x(n)

t → t + kh print k, t, x end for

19

Accuracy of Taylor method

• How many digits of the approximation are

reliable?

• Coarse assessment of the error for n = 4:
• The first term not included in Taylor series

is 1

5!h5x(5)(t)

• For example, if h =

1 100 the value of h5 is

1010 and depending on the order of mag- nitude of x(5) we can assess the error.

20

• Two types of errors are to be considered

here: – Local truncation error: Error we make when computing x(t+h) from x(t), this is estimated by

1 5!h5x(5) or, in general,

• rder O(hn)

– The second type of error is the cumu- lated effect of the iteration (x(t) is al- ready wrong, for t > a)

Roundoff error

• In addition the roundoff error has to be

considered

• In one iteration roundoff errors are typically

small, but accumulated over several itera- tions their impact can grow

• Interval arithmetics can be used to round-
• ff error

21

Interval arithmetics

• Some rules of interval arithmetics, e.g. for

a, b, c, d > 0: [a, b] + [c, d] = [a + c, b + d] (13) [a, b] + [c, d] = [a − d, b − c] (14) [a, b] ∗ [c, d] = [a ∗ c, b ∗ d] (15) [a, b]/[c, d] = [a/d, b/d] (16)

• For bounding the result of the error f([a, b])

the inclusion function can be derived with f([a, b]) = [min{f(x)|x ∈ [a, b]}, max{f(x)|x ∈ [a, b]}]

• Roundoff errors are particular harmful in

iterated systems with high Ljapunov coef- ficient

22

Runge Kutta method

• Runge Kutta methods∗ imitate the Taylor

method of order n without needing an ex- plicit formula for higher order derivatives

• Order 4 Runge Kutta method reads:

x(t + h) = x(t) + 1 6(F1 + 2F2 + 2F3 + F4) where F1 = hf(t, x) (17) F2 = hf(t + 1 2h, x + 1 2F1) (18) F3 = hf(t + 1 2h, x + 1 2F2) (19) F4 = hf(t + h, x + F3) (20)

• the next iterate is computed at the expense
• f computing f four times

∗Carl Runge and William Kutta

23

Runge Kutta Method (2)

• The error of the order-4 Runge Kutta method

is of order O(h4), as opposed to the order- 4 Taylor method with error of order O(h5)

• In practical situations a tolerance is pre-

scribed and the method must not produce an error beyond this tolerance

• The adaptive Runge Kutta Fehlberg method

adapts the stepsize h to reach the desired accuracy

• The allowable step-size is estimated in each

step

• Depending on whether the bound is met

the step-size is either increased or decreased

24

Simulating a system of ODEs X =

  

x1 . . . xm

   , F =   

f1 . . . fd

   , F ′ =   

f′

1

. . . f′

d

  

Now, we can approximate using Taylor poly- nominal: X(t + h) ≈ X(t) + hX′(t) + h21 2X′′(t) + ... Euler and Taylor methods are straightforward to apply. Runge-Kutta reads: X(t + h) = X + h 6(F1 + 2F2 + 2F3 + F4) F1 = F(X) (21) F2 = F(X + 1 2hF1) (22) F3 = F(X + 1 2hF2) (23) F4 = F(X + hF3) (24)

25

Higher order differential equations x(n) = f(t, x, x′, x′′, . . . , x(n−1)) with given initial value x(a), x′(a), . . . , xn−1(a) This DGL is rewritten as a system of coupled ODEs using: x0 = t, x1 = x, x2 = x′, x3 = x′′, . . . , xn = x(n−1) and x′ = 1 x′

1

= x2 x′

2

= x3 . . . x′

n−1

= xn x′

n

= f(x0, x1, . . . , xn) (25) Now, we can simulate ODE using systems sim- ulation.

26

Simulating Partial Differential Equations

• Different classical types of PDE are distin-

guished

• Often either a complete solution in a time

step, or boundary conditions that are pre- scribed for the perimeter are provided

• If there exists not a unique solution to the

PDE, we say the PDE is ill-posed

• finite difference methods and finite elements

methods are widely used to solve

• Finite difference methods for initial value

problems use local neighborhood to estab- lish solution of PDE in next time step

27

Example: Heat diffusion in 1-D ∂2u ∂x2 = ∂ ∂tu(x, t) u(0, t) = u(1, t) = 0 u(x, 0) = sin(πx) We can write f′(x) = 1 h(f(x + h) − f(x)) f′′(x) = 1 h2(f(x + h) − 2f(x) + f(x − h)) We rewrite the ODE using two different step length h (space),k (time): 1 h2(u(x + h, t) − 2u(x, t) + u(x − h, t)) = 1 k(u(x, t + k) − u(x, t)) We solve for u to obtain iteration rule: u(x, t + k) = k h2u(x + h, t) + (1 − 2 k h2)u(x, t) + k h2u(x − h, t)

28

Boundary value problem in 2-D Consider Laplacian equation: ∇2u = ∂2u ∂x2 + ∂2u ∂y2 = 0 with u(x, y) known at the boundary of the inte- gration region R. Five point formula (like ’von Neumann’ neighborhood): ∇2u ≈ 1 h2(u(x + h, y) + u(x − h, y) +u(x, y + h) + u(x, y − h)) − 4u(x, y) The inherent local error of the 5-point formula is in the order O(h2). Alternatively a 9 point formula can be used (moore neighborhood).

29

An large equation system on a grid can be used to solve the PDE: Consider the quadratic 5 × 5 grid xi = ih and yj = jh, i, j ∈ {1, . . . , 5}. Then the equation system reads: u12 + u32 + u23 + u21 − 4u22 = u22 + u42 + u33 + u31 − 4u32 = u52 + u32 + u43 + u41 − 4u42 = u13 + u33 + u24 + u22 − 4u23 = u43 + u23 + u34 + u32 − 4u33 = u53 + u33 + u44 + u42 − 4u43 = u34 + u14 + u25 + u23 − 4u24 = u44 + u24 + u35 + u33 − 4u34 = u54 + u34 + u45 + u43 − 4u44 = We use abbreviation uij = u(xi, yj). Gauss-Seidel method or sparse matrix tech- niques are employed to solve for 9 unknowns uij, with 1 < i, j < 5.

Hyperbolic problems (vibrating string)

• The given examples were a so-called parabolic

problem (heat diffusion) and an elliptic prob- lem (laplace equation).

• A third important type is the hyperbolic

type of equations, represented by the wave equation.

• In one spatial variable it reads:

∂2u ∂x2 = ∂2u ∂t2

• Typical boundary condition and initial value

are: u(x, 0) = f(x), ∂u ∂t = 0 u(0, t) = u(1, t) = 0

30

• For this problem, it can be derived analyt-

ically that: u(x, t) = 0.5(f(x + t) + f(x − t))

• provided f is twice differentiable and:

f(−x) = −f(x), f(x + 2) = f(x)

• this makes numerical computation very easy
• structure of problem is of great importance,

when solving PDE.

Summary DE

• Differential Equations can be used to model

continuous system dynamics

• They can be solved numerically and ana-

lytically (can be very difficult)

• Coupled differential equations introduce vec-

tor valued states

• Partial differential equations consider spatio-

temporal phenomena, such as diffusion and flow

• Initial values, and or boundary conditions,

need to be defined to obtain a unique so- lution (trajectory)

31

• Euler method is most simple method and

based on linear Taylor approximation

• Higher order Taylor approximation can pro-

vide more accurate methods

• Runge Kutta method consideres also the

numerical estimation of derivatives Literature: Ward Cheney and David Kincaid: Numerical Mathematics and Computing, 3rd edition, Brooks/Cole, CA, 1994

Summary Simulation Models

• Classification in state space, time discretiza-

tion, stochasticity

• Iterated function systems: continuous state,

discrete time

• Lindenmayer systems: string sequance state

space, discrete time

• Cellular Automaton: discrete cellular state

space, discrete time

• Markov Chains: discrete state, discrete time,

stochastic

• ODE: continuous state, continuous time

32

• Coupled ODE: continuous vector state space,

continuous time

• Partial DE: continuous space - spatio-temporal,

continuous time

• Stochastic DE: continuous space, contin-

uous time, stochastic (not covered)

• Discrete Event Simulation: continuous space,

discrete time (not covered)