Pdes
Pdes
Pdes
Hans Fangohr
May 3, 2012
1/
Outline I
3 Diffusion equation
6 Finite elements
7 Summary
2/
This lecture
3/
What are partial differential equations (PDEs)
Ordinary Differential Equations (ODEs)
one independent variable, for example t in
d 2x k
=− x
dt2 m
often the indepent variable t is the
time solution is function x(t)
important for dynamical systems, population growth,
control, moving particles
Partial Differential Equations (ODEs)
multiple independent variables, for example t, x and y in
∂u ∂2u ∂2u = +
D ∂x2 ∂y2
∂t
solution is function u(t, x, y)
important for fluid dynamics, chemistry,
electromagnetism, . . . , generally problems with spatial
resolution
4/
2d Diffusion equation
∂u ∂ 2u ∂ 2u = +
D ∂x2 ∂y2
∂t
u(t, x, y) is the concentration [mol/m3]
t is the time [s]
x is the x-coordinate [m]
y is the y-coordinate [m]
D is the diffusion coefficient [m2/s]
Also known as Fick’s second law. The heat equation has
the same structure (and u represents the temperature).
Example:
http://www.youtube.com/watch?v=WC6Kj5ySWkQ
5/
Examples of PDEs
Cahn Hilliard Equation (phase separation)
7/
Overview
Motivation:
We need derivatives of functions for example for
optimisation and root finding algorithms
Not always is the function analytically known (but we
are usually able to compute the function numerically)
The material presented here forms the basis of
the finite-difference technique that is commonly
used to solve ordinary and partial differential
equations.
The following slides show
the forward difference technique
the backward difference technique and the
central difference technique to approximate the
derivative of a function.
We also derive the accuracy of each of these methods.
8/
The 1st
derivative
(Possible) Definition of the derivative (or “differential
operator” dxd )
r df f (x + h) − f (x)
f (x) ≡ dx (x) = h→0
lim h
Use difference operator to approximate differential
operator
r df f (x + h) − f (x) f (x + h) − f
(x)
f (x) = (x) = lim ≈
dx h→0 h h
⇒ can now compute an approximation of f (x) simply by
r
evaluating f (twice).
This is called the forward difference because we use f (x)
and f (x + h).
9/
The 1st
derivative
Important questions: How accurate is this approximation?
10 /
Accuracy of the forward difference
Formal derivation using the Taylor series of f around x
f (x + h) ∞
Σ f (n)
(x)
= hn
n=0
n!
frr(x) f rrr(x)
= f (x) + hfr(x) + h2 + h3 +...
2! 3!
Rearranging for fr(x)
frr(x) f rrr(x)
hfr(x) = f (x + h) − f (x) − h2 − h3 − ...
2! 3!
1 frr(x) f rrr(x)
fr(x) = f (x + h) − f (x) − h 2
− h3 − ...
h 2! 3!
f (x + h) − f
(x) f′′(x) f′′′(x)
h2 2! − h3 3!
= − − ...
h h
11 / 61
Accuracy of the forward difference
f (x + h) − f (x) frr(x) 2 f (x)
rrr
= h − h 2! − h 3! − . . .
12 / 61
Accuracy of the forward difference (2)
r f (x + h) − f (x)
f (x) = + O(h)
13 / 61
Accuracy of the forward difference (2)
h
14 / 61
The 1st derivative using the backward difference
16 / 61
Accuracy of the backward difference
Formal derivation using the Taylor Series of f around x
frr(x) f rrr(x)
f (x − h) = f (x) − hf (x) + h
r 2
−h 3
+...
2! 3!
Rearranging for f (x)
r
frr(x) f rrr(x)
hfr(x) = f (x) − f (x − h) + h2 − h3 − ...
2! 3!
1 frr(x) f rrr(x)
f (x)
r
= f (x) − f (x − h) + h 2
−h 3
− ...
h 2! 3!
f (x) − f (x −
h) f′′(x)
h2 2! − h3 3!
f′′′(x)
= + − ...
h h
f (x) − f (x − h) frr(x) 2 f (x)
rrr
= +h −h − ...
17 / 61
Accuracy of the backward difference
h 2! 3!
18 / 61
Accuracy of the backward difference (2)
f (x) − f (x − h)
r
f (x) = + back (h) (1)
E
h
r f (x) − f (x − h)
f (x) = + O(h) 19 / 61
Accuracy of the backward difference (2)
h
20 / 61
Combining backward and forward differences (1)
The approximations are
forward:
f (x + h) − f (x) (2)
+ forw (h)
r
f (x) =
E
backward h
(h) (3)
r f (x) − f (x − h)
f (x) =
E +
back
h
frr(x) f rrr(x) f rrrr(x) f rrrrr(x)
E (h) = −h −h 2
−h 3
−h 4
− ...
forw
2! 3! 4! 5!
frr(x) 2 f (x)
rrr
3f
rrrr
(x) 4f
rrrrr
(x)
E back (h) = h − h +h − h +...
2! 3! 4! 5!
21 / 61
Combining backward and forward differences (1)
⇒ Add equations (2) and (3) together, then the error
cancels partly!
22 / 61
Combining backward and forward differences (2)
Add these lines together
r f (x + h) − f (x)
f (x) = + forw (h)
E
h (h)
r f (x) − f (x − h)
f (x) =
E +
back
h
f (x + h) − f (x − h)
+ forw(h) + E back (h)
r
2f (x) =
E
h
Adding the error terms:
E (h) + E f rrr(x) f rrrrr(x)
(h) = −2h 2
− 2h 4
− ...
forw back
3! 5!
The combined (central) difference operator is h)
fh +
Combining backward and forward differences (2)
E 2h
cent (h)
with f rrr(x) f rrrrr(x)
E (h) = −h2 − h4 − ...
cent
3! 5! 16 / 47
Central difference
Can be derived (as on previous slides) by adding forward
and backward difference
Can also be interpreted geometrically by defining the
differential operator as
df
(x) = lim f (x + h) − f (x −
dx h) 2h
h→0
18 / 61
Example (1)
20 / 61
Example (2)
80
60
40
20
0
0 1 2 3 4 5
x
21 / 61
Summary
22 / 61
Appendix: source to compute figure on page 19
I
EPS =1 # very large EPS to provoke inaccuracy
23 / 61
Appendix: source to compute figure on page 19
II def f( x):
""" Our test funtion with
convenient property that
df/ dx = f"""
return np. exp ( x)
24 / 61
Appendix: source to compute figure on page 19
III
25 / 61
Note: Euler’s (integration) method — derivation
using finite difference operator
Use forward difference operator to approximate
differential operator
dy y(x + h) − y(x)
(x) = lim y(x + h) −
dx ≈ y(x) h
h→0 h
Change differential to difference operator in dy
= f (x, y)
dx
dy
f (x, y) = y(x + h) −
y(x) h
≈
dx
hf (x, y) ≈ y(x + h) − y(x)
=⇒ yi+1 = yi + hf (xi, yi)
26 / 61
⇒ Euler’s method (for ODEs) can be derived from the
forward difference operator.
27 / 61
Note: Newton’s (root finding) method —
derivation from Taylor series
We are looking for a root, i.e. we are looking for a x so that
f (x) = 0.
We have an initial guess x0 which we refine in subsequent iterations:
f (xi )
xi+1 = — where hi = . (4)
. xi hi f′(x i)
29 / 61
Diffusion equation
2
The 2d operator ∂
∂22 is called the Laplace
+ ∂y
∂x2
operator ∆, so that we can also write
∂u ∂2 + 2
= ∂2
D ∂x2 ∂y = D∆u
∂t
The diffusion equation (with constant diffusion
coefficient D) reads ∂t∂u = D∆u where the Laplace
operator depends on the number d of spatial dimensions
2
∂
d = 1: ∆ = ∂x2
2
d = 2: ∆ = ∂
+ ∂2
∂x2 ∂y2
2
d = 3: ∆ = ∂
∂x2 + ∂2 + ∂2
∂y2 ∂z2
30 / 61
2
1d Diffusion equation∂t∂u = D∂∂xu2
In one spatial dimension, the diffusion equation reads
∂u ∂ 2u
=D
∂t ∂x2
This is the equation we will use as an example.
Let’s assume an initial concentration
u(x, t0) = √1 exp −(x—xmean )2
σ σ 2 with x mean = 0 and
2π
width σ = 0.5.
0.8
0.7 (x,t0 )
0.6
0.5
u(x,t_0)
0.4
0.3
0.2
0.1
0.0
6 4 2 0 2 4 6
x 31 / 61
2
∂u ∂ u
1d Diffusion eqn ∂t
= D ∂x2, time integration
I
Let us assume that we have some way of computing
2
∂u
D ∂x 2 at time and let’s call this g(x, t0), i.e.
t0
∂2u(x, t0)
g(x, t0) ≡ D
∂x2
We like to solve
∂u(x, t)
= g(x, t0)
∂t
to compute u(x, t1) at some later time t1.
Use finite difference time integration scheme:
Introduce a time step size h so that t1 = t0 + h.
32 / 61
2u
∂u
1d Diffusion eqn ∂t
= D ∂x2, time integration
∂
II
Change differential operator to forward difference operator
∂u(x, t)
g(x, t0) = u(x, t0 + h) − u(x, t0)
∂t = lim
h→0 h (6)
u(x, t0 + h) − u(x, t0) (7)
≈ h
Rearrange to find u(x, t1) ≡ u(x, t0 + h) gives
I
∂u ∂ 2u
=D = g(x, t)
∂t ∂x 2
2
u(x,t)
Need to compute g(x, t) = D∂ for a given u(x, t).
∂x2
Can ignore the time dependence here, and obtain
∂2u(x)
g(x) = D .
∂x2
Recall
that ∂ 2u ∂ ∂u
2
= ∂u
∂x ∂x ∂x
and we that know how to compute∂x
34 / 61
2u
∂u
1d Diffusion eqn ∂t
= D ∂x2, spatial part
∂
I using central
differences.
35 / 61
Second order derivatives from finite differences I
dx h
36 / 61
Second order derivatives from finite differences II
2
f
Apply the central difference equation twice to obtaindxd2 :
d2f d df
2
(x) = dx dx(x)
dx
2 h) − f (x −
d f (x + 1 2 1 h)
≈
dx h
1 d 1 d 1
= f x+ h − f x− h
h dx 2 dx 2
1 f ( x + h) − f ( x ) f ( x ) − f ( x − h )
≈ h h − h
f (x + h) − 2f (x) + f (x − h)
= (9)
h2
37 / 61
2
Recipe to solve∂t∂u = D∂∂xu2
38 / 61
2u
∂u
A sample solution ∂t
= D ∂x2, I
∂
import numpy as np
import matplotlib . pyplot as plt
import matplotlib . animation as animation
39 / 61
2u
∂u
A sample solution ∂t
= D ∂x2, II
∂
def compute_g ( u, D, h ):
""" given a u(x, t) in array , compute g(x, t)=D* d^2 u/ dx^2
using central differences with spacing h,
and return g(x, t). """
d 2 u_dx2 = np. zeros( u. shape , np. float)
for i in range (1 , len ( u ) -1):
d 2 u_dx2 [ i] = ( u[ i+1] - 2* u[ i]+u[i -1])/ h
**2 # special cases at boundary: assume Neuman
boundary # conditions , i. e. no change of u over
boundary
# so that u[0]- u [ - 1]=0 and thus u [ - 1]=u[0]
i=0
d 2 u_dx2 [ i] = ( u[ i+1] - 2* u[ i]+u[ i])/ h
**2 # same at other end so that u[N -1]- u[
N]=0 # and thus u[ N]=u[N -1]
i= len ( u)- 1
d 2 u_dx2 [ i] = ( u[ i] - 2* u[ i]+u[i -1])/ h
**2 return D* d 2 u_dx2
40 / 61
2u
∂u
A sample solution ∂t
= D ∂x2, II
∂
u = u +dt* g
41 / 61
2u
∂u
A sample solution ∂t
= D ∂x2, III
∂
return u
42 / 61
2u
∂u
A sample solution ∂t
= D ∂x2, III
∂
43 / 61
2u
∂u
A sample solution ∂t
= D ∂x2, IV
∂
return l,
44 / 61
Boundary conditions I
46 / 61
Numerical issues
47 / 61
Performance issues
48 / 61
Finite Elements
Another widely spread way of solving PDEs is using so-called
finite elements.
Mathematically, the solution u(x) for a problem like
∂ 2u
∂x2
= f (x) is written as
N
Σ
u(x) = uiφi(x) (10)
i=1
where each ui is a number (a coefficient), and each φi(x)
a known function of space.
The φi are called basis or shape functions.
Each φi is normally chosen to be zero for nearly all x,
and to be non-zero close to a particular node in the
finite element mesh.
By substitution (10) into the PDE, a matrix system
can be obtained, which – if solved – provides the
coefficients ui, and thus the solution.
49 / 61
Finite Elements vs Finite differences
Finite differences
are mathematically much simpler and
for simple geometries (such as cuboids) easier to
program
Finite elements
have greater flexibility in the shape of the
domain, the specification and implementation of
boundary
conditions is easier
but the basic mathematics and code is more
complicated.
50 / 61
Practical observation on time integration
51 / 61
Summary
52 / 61
changeset: 53:a22b7f13329e
tag: tip
user: Hans Fangohr [MBP13] <fangohr@soton.ac.uk>
date: Fri Dec 16 10:57:15 2011 +0000
53 / 61