Introduction to Fourier Transforms for Solving Partial Differential

Introduction to Fourier Transforms for Solving
Partial Differential Equations and a
Supplement on Monte Carlo Methods
Benson Muite
[email protected]
http://math.ut.ee/˜benson
http://en.wikibooks.org/wiki/Parallel_Spectral_Numerical_Methods
26 April 2014
Outline
• Fourier Series
• The Heat Equation
• The Allen Cahn Equation
• Splitting
• Iteration
• Parallel Computing
• Monte Carlo Methods for Computing Integrals
Fourier Series: Separation of Variables 1
dy
dt
dy
y
Z
dy
y
ln y + a
e
ln y+a
= y
= dt
Z
=
dt
= t +b
= et+b
eln y ea = et eb
eb t
e
y =
ea
y (t) = cet
Fourier Series: Separation of Variables 2
•
ut = uxx
• Suppose u = X (x)T (t)
•
dT
dt (t)
T (t)
=
d2 X
(x)
dx 2
X (x)
= −C,
• Solving each of these separately and then using linearity
we get a general solution
•
∞
X
n=0
αn exp(−Cn t) sin(
p
p
Cn x) + βn exp(−Cn t) cos( Cn x)
Fourier Series: Separation of Variables 3
• How do we find a particular solution?
• Suppose u(x, t = 0) = f (x)
• Suppose u(0, t) = u(2π, t) and ux (0, t) = ux (2π, t) then
recall
•
Z
2π
π m=n
,
0 m=
6 n
0
Z 2π
π m=n
,
cos(nx) cos(mx) =
0 m 6= n
0
Z 2π
cos(nx) sin(mx) = 0.
sin(nx) sin(mx) =
0
Fourier Series: Separation of Variables 4
• So if
f (x) =
∞
X
αn sin(nx) + βn cos(nx).
n=0
• then
R 2π
αn =
f (x) sin(nx)dx
R 2π
2
0 sin (nx)dx
0
R 2π
βn =
f (x) cos(nx)dx
.
R 2π
2
0 cos (nx)dx
0
• and
u(x, t) =
∞
X
n=0
exp(−n2 t) [αn sin(nx) + βn cos(nx)]
Fourier Series: Separation of Variables 5
• The Fast Fourier Transform allows one to find good
approximations to αn and βn when the solution is found at
a finite number of evenly spaced grid points
• By rescaling, can consider intervals other than [0, 2π)
• Fourier transform also works on infinite intervals, but
require function to decay to the same constant value at ±∞
Complex Fourier Series
• By using Euler’s formula, one can get a simpler expression
for a Fourier series where sine and cosine are combined
P
• u= ∞
n=−∞ γn exp(inx) x ∈ [0, 2π)
Source: http://en.wikipedia.org/wiki/File:Euler%27s_formula.svg
A Computational Algorithm for Computing An
Approximate Fourier Transform 1
• Analytic method of computing Fourier transform can be
tedious
• Can use quadrature to numerically evaluate Fourier
transforms – O(n2 ) operations
• Gauss and then Cooley and Tukey found O(n log n)
algorithm
• Key observation is to use factorization and recursion
• Modern computers use variants of this idea that are more
suitable for computer hardware where moving data is more
expensive than floating point operations
A Computational Algorithm for Computing An
Approximate Fourier Transform 2
Example pseudo code to compute a radix 2 out of place DFT where x has length that is a power of 2
1: procedure X0,...,N−1 (ditfft2(x,N,s))
2: DFT of (x0 , xs , x2 s, ..., x(N−1)s )
3: if N=1 then
4:
trivial size-1 DFT of base case
5:
X0 ← X0
6: else
7:
DFT of (x0 , x2s , x4s , ...)
8:
X0,...,N/2−1 ← ditfft2(x, N/2, 2s)
9:
DFT of (xs , xs+2s , xs+4s , ...)
10:
XN/2,...,N−1 ← ditfft2(x + s, N/2, 2s)
11:
Combine DFTs of two halves into full DFT
12:
for k = 0 → N/2 − 1 do
13:
t ← Xk
14:
Xk ← t + exp(−2πik /N)Xk +N/2
15:
Xk +N/2 ← t − exp(−2πik/N)Xk +N/2
16:
end for
17: end if
18: end procedure
Sources: http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm,
http://cnx.org/content/m16336/latest/
The Heat Equation: Finding Derivatives and
Timestepping
• Let
u(x) =
X
ˆk exp(ikx)
u
k
• then
dν u X
ˆk .
=
(ik)ν u
dx ν
• Consider ut = uxx , which is approximated by
ˆk
∂u
ˆk
= (ik)2 u
∂t
ˆkn+1 − u
ˆkn
u
ˆkn+1
= (ik)2 u
δt
ˆkn
ˆkn+1 [1 − δt(ik)2 ] = u
u
ˆkn
u
ˆkn+1 =
.
u
[1 − δt(ik )2 ]
The Heat Equation: Finding Derivatives and
Timestepping
• Python demonstration
The Allen Cahn Equation: Implicit-Explicit Method
• Consider ut = uxx + u − u 3 , which is approximated by
backward Euler for the linear terms and forward Euler for
the nonlinear terms
ˆ
∂u
∂t
ˆ n+1 − u
ˆn
u
δt
ˆ+u
ˆ − uˆ3
= (ik)2 u
ˆ + uˆn − (uˆn )3
= (ik)2 u n+1
The Allen Cahn Equation: Implicit-Explicit Method
• Python demonstration
The Allen Cahn Equation: Splitting
• Again consider ut = uxx + u − u 3
• Can solve the equation easily using Fourier Series
ˆ = exp
ut = uxx + u, u
ˆ0
(ik)2 + 1 t u
• Can solve the ordinary differential equation ut = u 3 exactly,
u = √ u2o
uo t+1
• Idea of splitting is to solve the two equations exactly in
alternating steps
• Solve ut = uxx + u for a time δt, take that solution as an
initial condition and then solve ut = u 3 for a time δt.
• Note that the splitting ut = uxx and ut = u − u 3 is also
possible, but care is needed in solving the nonlinear
equation
The Allen Cahn Equation: Splitting
Example pseudo code to explain Godunov splitting
1: procedure G ODUNOV A LLEN C AHN(u)
2:
u = uo
3:
for n = 0 → Nt do
ˆ ← FFT (u)
4:
u
ˆ ← exp([(ik)2 + 1]δt)u
ˆ
5:
u
ˆ
6:
u ← IFFT (u )
7:
u ← √ 2u
u δt+1
8:
end for
9: end procedure
The Allen Cahn Equation: Splitting
• Python Demonstration
The Viscous Burgers Equation: Iteration
• Consider using splitting for ut = 0.5(u 2 )x + uxx
• Solve exactly ut = uxx using Fourier series
• Solve ut = uux using backward Euler
2 u n+1,k+1 − u n
n+1,k
= 0.5 u
δt
x
• Choose u n+1,1 = u n and stop iterating on k when
|u n+1,k+1 − u n+1,k |
or
n+1,k+1
2 u
− un
n+1,k+1
− u
δt
x
is small enough
Reference: http://en.wikipedia.org/wiki/Burgers%27_equation
The Viscous Burgers Equation: Iteration
• Python Demonstration
Parallel Computing
• Due to limitations in device physics
power = resistance × current2
to get faster results, instead of increasing the frequency
and hence the current, want to have many devices at same
current, so get a linear increase in power, keeping all other
things fixed
• For more on this, take distributed computing course next
semester. We give a brief glimpse here.
• Will use Rocket, Fortran (primarily Fortran90 though
knowledge of Fortran77 is still useful) and MPI, though
similar ideas and methods work on other parallel
computers such as your cell phone, laptop, video game
console etc.
Parallel Computing
• Main idea is to split up work into independent tasks with as
few synchronisations as possible
• Kalev, Kamau and Khruschev (or Svetlana, Tiia and
Wanjiku) need to make 99 salads
• Possible organizational options:
• Kalev, Kamau and Khruschev each make 33 salads.
• Simultaneously Kalev chops beetroots, Kamau chops
tomatoes, Khruschev chops carrots. They then mix the
ingredients in one bowl.
• Can you come up with others? What is the fastest
assuming each of the three people are identically good at
all the tasks? What if some are better at some tasks than
others?
Parallel Computing: MPI
• MPI (Message Passing Interface) - one way of making a
program parallel
• Standard can be found at:
http://www.mpi-forum.org
• Use a library to allow computers to talk to each other by
sending messages and having some explicit co-ordination
• MPI works for C and Fortran programs. Some unofficial
bindings available for other programs, such as Python:
http://mpi4py.scipy.org/
• Many online resources available including:
https://computing.llnl.gov/tutorials/mpi/
http://www.shodor.org/refdesk/Resources/
Tutorials/BasicMPI/
Hello World Example Program
https://github.com/openmichigan/PSNM/blob/
master/IntroductionToParallelProgramming/
Programs/HelloworldMpi/helloworld.f90
Hello World Example Rocket Submission Script
# ! / b i n / bash
#SBATCH −N 2
#SBATCH −−sockets−per−node=2
#SBATCH −−cores−per−s o c k e t =10
#SBATCH −−threads−per−core =1
#SBATCH −t 0 0 : 1 0: 0 0
module l o a d i n t e l c l u s t e r s t u d i o
export I MPI PMI LIBRARY = / u s r / l i b 6 4 / l i b p m i . so
srun −−c p u b i n d =cores
helloworld
Listing 1: An example submission script for use on Rocket.
Monte Carlo Method: A Probabilistic Way to Calculate
Integrals
• Recall
¯f =
1
b−a
Z
b
f (x) dx
a
• Hence given ¯
f
Z
b
f (x)dx = (b − a)¯f
a
• Doing the same in 2 dimensions and estimating the error
using the standard deviation
s
ZZ
f (x, y ) dA ≈ A(R)¯f ± A(R)
f 2 − (¯f )2
,
N
R
• Approximate ¯
f by random sampling
¯f ≈
PN
i=1 f (xi , yi )
N
PN
and
f2
≈
2
i=1 (f (xi , yi ))
N
Monte Carlo Method: Python Program
”””
A program t o approximate an i n t e g r a l u s i n g a Monte C a r l o method
T h i s c o u l d be made f a s t e r by u s i n g v e c t o r i z a t i o n , however i t i s
k e p t as s i m p l e as p o s s i b l e f o r c l a r i t y and ease o f t r a n s l a t i o n i n t o
o t h e r languages
”””
import math
import numpy
import t i m e
numpoints =4096 # number o f random sample p o i n t s
I 2 d =0.0
# i n i t i a l i z e value
I2dsquare =0.0
# i n i t i a l i z e to allow f o r c a l c u l a t i o n of variance
f o r n i n xrange ( numpoints ) :
x=numpy . random . u n i f o r m ( )
y =4.0∗numpy . random . u n i f o r m ( )
I 2 d = I 2 d +x∗x +2.0∗ y∗y
I2dsquare = I2dsquare +( x∗x +2.0∗ y∗y)∗∗2
# we s c a l e t h e i n t e g r a l by t h e t o t a l area and d i v i d e by t h e number o f
# p o i n t s used
I 2 d = I 2 d ∗4/ numpoints
I2dsquare = I2dsquare ∗4/ numpoints
E s t i m E r r o r =4∗numpy . s q r t ( ( I 2 d∗∗2−I2dsquare ) / numpoints ) # e s t i m a t e d e r r o r
p r i n t ” Value : %f ” %I 2 d
p r i n t ” E r r o r e s t i m a t e : %f ” %E s t i m E r r o r
Listing 2: A Python program which demonstrates how to use the
Monte Carlo method to calculate the volume below z = x 2 +2y 2 ,
with (x, y) ∈ (0, 1) × (0, 4).
Sample Results of Monte Carlo Program
N
16
256
4096
65536
∞
Value
41.3026
47.1855
43.4527
44.0026
44
Error Estimate
+/- 30.9791
+/- 9.0386
+/- 2.0280
+/- 0.5151
0.0
Monte Carlo Method: Serial Fortran Program
https://github.com/openmichigan/PSNM/blob/
master/IntroductionToParallelProgramming/
Programs/montecarloserial/montecarloserial.f90
Monte Carlo Method: Parallel Fortran Program
https:
//github.com/openmichigan/PSNM/blob/master/
IntroductionToParallelProgramming/Programs/
montecarloparallel/montecarloparallel.f90
New Key Concepts
• Discrete Fourier Transform
• Fast Fourier Transform
• Splitting Methods – take a hard problem and make it into
smaller easier problems
• Iteration – when you cannot solve exactly, solve
approximately to a good tolerance
• Parallel Computing – MPI
• Monte Carlo Method – method of calculating integrals
Further Work
• Several projects possible, either as summer projects or
coursework
• Some of these are extensions of the homework, consider
trying the optional exercises
• Please also contact me if interested in building a
SUPERCOMPUTER
References
• Cooley J.W. and Tukey J.W. “An algorithm for the machine calculation of complex
Fourier series” Math. Comput. 19, 297–301, (1965)
• Heideman M.T., Johnson D.H. and Burrus C.S. “Gauss and the History of the
Fast Fourier Transform” IEEE ASSP Magazine 1, 14–21, (1984)
¨
• Vainikko E. Fortran 95 Ja MPI Tartu Ulikool
Kirjastus (2004)
http://kodu.ut.ee/˜eero/PC/F95jaMPI.pdf
• Chen G., Cloutier B., Li N., Muite B.K., Rigge P. and Balakrishnan S., Souza A.,
West J. “Parallel Spectral Numerical Methods” http://shodor.org/
petascale/materials/UPModules/Parallel_Spectral_Methods/
• Corral M. Vector Calculus http://www.mecmath.net/
• Greenbaum A. and Chartier T.P. Numerical Methods Princeton University Press
(2012)
• Sauer T. Numerical Analysis Pearson (2012)
• Petersen W.P. and Arbenz P. Introduction to Parallel Computing Oxford University
Press (2004)
Acknowledgements
• Oleg Batrashev
• Leonid Dorogin
• Michael Quell
• Eero Vainikko
• The Blue Waters Undergraduate Petascale Education Program
administered by the Shodor foundation
• The Division of Literature, Sciences and Arts at the University of
Michigan
Fun with Fourier Series
• A Fourier series can represent a regular k-gon in the
complex plane
P
−2 exp [i(1 + kn)t]
• f (t) = +∞
n=−∞ (1 + kn)
t ∈ (−π, π)
Robert, A “Fourier Series of Polygons” Amer. Math. Monthly 101(5) pp. 420-428 (1994)
Schonberg, IJ “The Finite Fourier Series And Elementary Geometry” Amer. Math. Monthly 57(6) pp. 390-404 (1950)

Download Report