Pdes

Download as docx, pdf, or txt
Download as docx, pdf, or txt
You are on page 1of 61

Solving partial differential equations (PDEs)

Hans Fangohr

Engineering and the Environment


University of Southampton
United Kingdom
fangohr@soton.ac.uk

May 3, 2012

1/
Outline I

1 Introduction: what are PDEs?

2 Computing derivatives using finite differences

3 Diffusion equation

4 Recipe to solve 1d diffusion equation

5 Boundary conditions, numerics, performance

6 Finite elements

7 Summary

2/
This lecture

tries to compress several years of material into 45 minutes


has lecture notes and code available for download at
http://www.soton.ac.uk/∼fangohr/teaching/comp6024

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)

Fluid dynamics (including ocean and atmospheric models,


plasma physics, gas turbine and aircraft modelling)
Structural mechanics and vibrations, superconductivity,
micromagnetics, . . .
6/
Computing derivatives
using finite differences

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) frr(x) 2 f (x)


rrr
f (x) = h − h 2! − h 3! − . . .
`
˛ ¸ x
Efo rw (h)
r f (x + h) − f (x)
f (x) = +
E forw (h)
h

Therefore, the error term Eforw(h) is


frr(x) f rrr(x)
forw
E (h) = −h − h2 − ...
2! 3!
Can also be expressed as

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

Another definition of the derivative (or “differential


operator” dxd )
df
(x) = lim f (x) − f (x −
dx h) h
h→0

Use difference operator to approximate differential


operator
df
(x) = lim f (x) − f (x − f (x) − f (x − h)
dx h) h ≈ h
h→0

This is called the backward difference because we use


15 / 61
The 1st derivative using the backward difference
f (x) and f (x − h).
How accurate is 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)

r f (x) − f (x − h) frr(x) 2 f (x)


rrr
f (x) = h + h` 2! − h ˛ ¸ 3! − . x
..
Eba c k (h)

f (x) − f (x − h)
r
f (x) = + back (h) (1)
E
h

Therefore, the error term Eback(h) is


frr(x) f rrr(x)
back
E (h) = h −h 2
− ...
2! 3!
Can also be expressed as

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

and taking the finite difference form


df
(x) f (x + h) − f (x − h)
≈ 2h
dx
Error of the central difference is only O(h2), i.e. better
than forward or backward difference
17 / 61
Central difference

It is generally the case that symmetric differences


are more accurate than asymmetric
expressions.

18 / 61
Example (1)

Using forward difference to estimate the derivative of


f (x) = exp(x)

fr(x) ≈ fr f (x + h) − f (x) exp(x + h) − exp(x)


= =
forw
h h
Numerical example:
h = 0.1, x = 1
fr(1) ≈ (1.0) =
exp(1.1)—exp(1)
= 2.8588
r
ff orw 0.1
Exact answers is fr(1.0) = exp(1) = 2.71828
(Central diff: f′ 0.2
cent
19 / 61
Example (1)
(1.0) =
exp(1+0.1)−exp(1−0.1)
=
2.72281)

20 / 61
Example (2)

Comparison: forward difference, central difference and exact


derivative of f (x) = exp(x)
Approximations of df/dx for f(x)=exp(x)
140 forward h=1
120 central h=1
forward h=0.0001
100 exact
df/dx(x)

80

60

40

20

0
0 1 2 3 4 5
x

21 / 61
Summary

Can approximate derivatives of f numerically using only


function evaluations of f
size of step h very important
central differences has smallest error term

name formula error


f (x+h)—f (x)
forward fr(x) = h O(h)
backward fr(x) = f(x)—f(x—h)
h O(h)
f (x+h)—f (x— O(h2)
central
h ) fr(x) =
2h

22 / 61
Appendix: source to compute figure on page 19
I
EPS =1 # very large EPS to provoke inaccuracy

def forwarddiff (f,x, h= EPS ):


# df/ dx = ( f( x+h)- f( x) )/ h + O( h)
return ( f( x+h)- f( x) )/ h

def backwarddiff (f,x, h= EPS ):


# df/ dx = ( f( x)- f(x- h) )/ h + O( h)
return ( f( x)- f(x- h) )/ h

def centraldiff (f,x, h= EPS ):


# df/ dx = ( f( x+h) - f(x- h))/ h + O( h^2)
return ( f( x+h) - f(x- h ))/(2* h)

if name == " main ":


# create example plot
import pylab
import numpy as np
a=0 # left and
b=5 # right limits for x
N =11 # steps

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)

xs=np. linspace (a,b, N)


forward = []
forward_small_h = []
central = []
for x in xs:
forward . append ( forwarddiff (f, x)
) central. append ( centraldiff (f,
x) ) forward_small_h . append (
forwarddiff (f,x, h=1 e -4))

pylab . figure ( figsize =(6 ,4))


pylab . axis ([ a,b ,0 , np. exp ( b )])
pylab . plot( xs , forward ,’^’, label=’ forward h=% g’%
EPS) pylab . plot( xs , central ,’x’, label=’ central h=%
g’% EPS) pylab . plot( xs , forward_small_h ,’o’,
label=’ forward h=% g’% 1e -4)
xsfine = np. linspace (a,b, N *100)

24 / 61
Appendix: source to compute figure on page 19
III

pylab . plot( xsfine , f( xsfine ) , ’ - ’, label=’ exact ’)


pylab . grid ()
pylab . legend ( loc=’ upper left
’) pylab . xlabel(" x")
pylab . ylabel(" df/ dx( x)")
pylab . title (" Approximations of df/ dx for f( x)= exp (
x)") pylab . plot ()
pylab . savefig (’ central - and - forward - difference . pdf ’)
pylab . show ()

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)

This equation can be derived from the Taylor series of f around


x. Suppose we guess the root to be at x and x + h is the actual
location of the root (so h is unknown and f (x + h) = 0):
f (x + h) = f (x) + hf′(x) + . . .
0 = f (x) + hf′(x) + . . .
=⇒ 0 ≈ f (x) + hf′(x)
f (x)
⇐⇒ h ≈ − . (5)
f′(x)
28 / 61
The diffusion equation

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

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

u(x, t1) ≈ u(x, t0) + hg(x, t0)

We can generalise this using ti = t0 + ih to read

u(x, ti+1) ≈ u(x, ti) + hg(x, ti) (8)

→ If we can find g(x, ti), we can compute u(x, ti+1)


33 / 61
2u
∂u
1d Diffusion eqn ∂t
= D ∂x2, spatial part

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

Recall central difference equation for first order derivative


df
(x) f (x + h) − f (x − h)
≈ 2h
dx
will be more convenient to replace h by 1 h:
2
df f (x + 1 h) − f (x − 12 h)
(x) ≈ 2

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

1 Discretise solution u(x, t) into discrete values


2 uji ≡ u(xj, ti) where
xj ≡ x0 + j∆x and
ti ≡ t0 + i∆t.
3 Start with time iteration i = 0
4 Need to know configuration2 u(x, ti).
5 Then compute g(x, t ) = D ∂ u using finite differences
i 2
∂x
(9).
6 Then compute u(x, ti+1) based on g(x, ti) using (8)
7 increase i to i + 1, then go back to 5.

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

a, b=-5 ,5 # size of box


N = 51 # number of subdivisions
x=np. linspace (a,b, N) # positions of subdivisions
h=x[1]- x[0] # discretisation stepsize in x- direction

def total( u):


""" Computes total number of moles in u."""
return (( b- a)/ float( N)* np. sum ( u))

def gaussdistr ( mean , sigma , x):


""" Return gauss distribution for given numpy array x"""
return 1./( sigma * np. sqrt (2* np. pi ))* np. exp (
-0.5*( x- mean )**2/ sigma **2)

# starting configuration for u(x, t0 )


u = gaussdistr ( mean =0., sigma =0.5 , x=x)

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

def advance_time ( u, g, dt):


""" Given the array u, the rate of change array g,
and a timestep dt, compute the solution for u
after t, using simple Euler method ."""

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

# show example , quick and dirtly , lots of global


variables dt = 0.01 # step size or time
stepsbeforeupdatinggraph = 20 # plotting is slow
D = 1. # Diffusion coefficient
stepsdone = 0 # keep track of iterations

def do_steps(j, nsteps= stepsbeforeupdatinggraph ):


""" Function called by Func Animation class. Computes
nsteps iterations , i. e. carries forward solution
from u(x, t_i) to u(x, t_{ i+ nsteps }).
"""
global u, stepsdone
for i in range ( nsteps ):
g = compute_g ( u, D, h)
u = advance_time ( u, g, dt)
stepsdone += 1
time_passed = stepsdone * dt
print(" stepsdone =%5 d, time =%8 gs , total( u )=%8 g"
%
( stepsdone , time_passed , total( u )))
l. set_ydata ( u) # update data in plot

42 / 61
2u
∂u
A sample solution ∂t
= D ∂x2, III

fig1 . canvas. draw () # redraw the


canvas

43 / 61
2u
∂u
A sample solution ∂t
= D ∂x2, IV

return l,

fig1 = plt. figure () # setup animation


l,= plt. plot(x,u,’b- o’) # plot initial u(x, t)
# then compute solution and animate
line_ani = animation . FuncAnimation ( fig1 ,
do_steps , range (10000 ))
plt. show ()

44 / 61
Boundary conditions I

For ordinary differential equations (ODEs), we need to


know the initial value(s) to be able to compute a
solution.
For partial differential equations (PDEs), we need to
know the initial values and extra information about
the behaviour of the solution u(x, t) at the boundary
of the spatial domain (i.e. at x = a and x = b in this
example).
Commonly used boundary conditions are
Dirichlet boundary conditions: fix u(a) = c to some
constant.
Would correspond here to some mechanism that keeps
the concentration u at position x = a constant.
45 / 61
Boundary conditions II

Neuman boundary conditions: fix the change of u across


the boundary, i.e.
∂u
(a) = c.
∂x
For positive/negative c this corresponds to an imposed
concentration gradient.
For c = 0, this corresponds to conservation of the
atoms in the solution: as the gradient across the
boundary cannot change, no atoms can move out of
the box. (Used in our program on slide 35)

46 / 61
Numerical issues

The time integration scheme we use is explicit because we


have an explicit equation that tells us how to
compute u(x, ti+1) based on u(x, ti) (equation (8) on
slide 30)
An implicit scheme would compute u(x, ti+1) based on
u(x, ti) and on u(x, ti+1).
The implicit scheme is more complicated as it requires
solving an additional equation system just to find
u(x, ti+1) but allows larger step sizes ∆t for the time.
The explicit integration scheme becomes quickly unstable
if ∆t is too large. ∆t depends on the chose spatial
discretisation ∆x.

47 / 61
Performance issues

Our sample code is (nearly) as slow as possible


interpreted language
explicit for loops
enforced small step size from explicit scheme
Solutions:
Refactor for-loops into matrix operations and use
(compiled) matrix library (numpy for small systems, use
scipy.sparse for larger systems)
Use library function to carry out time integration (will
use implicit method if required), for example
scipy.integrate.odeint.

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

Usually, we solve the spatial part of a PDE using some


discretisation scheme such as finite differences and finite
elements).
This results in a set of coupled ordinary differential
equations (where time is the independent variable).
Can think of this as one ODE for every cube from
our discretisation.
This temporal part is then solved using time integration
schemes for (systems of) ordinary differential equations.

51 / 61
Summary

Partial differential equations important in many contexts


If no analytical solution known, use numerics.
Discretise the problem through
finite differences (replace differential with
difference operator, corresponds to chopping space
and time in little cuboids)
finite elements (project solution on localised basis
functions, often used with tetrahedral meshes)
related methods (finite volumes, meshless methods).
Finite elements and finite difference calculations are
standard tools in many areas of engineering, physics,
chemistry, but increasingly in other fields.

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

You might also like