Introduction to Numerical Methods for Solving Partial Differential Equations Benson Muite [email protected] http://kodu.ut.ee/˜benson http://en.wikibooks.org/wiki/Parallel_Spectral_Numerical_Methods 25 March 2015 Outline • Review • Aim • The heat equation • Methods for solving the heat equation • Fourier Series • The Allen Cahn Equation • Summary • References Ordinary Differential Equations: Numerical Schemes • Forward Euler method y n+1 − y n = f yn δt • Backward Euler method y n+1 − y n = f y n+1 δt • Implicit Midpoint rule n+1 y + yn y n+1 − y n =f δt 2 • Crank Nicolson Method f y n+1 + f (y n ) y n+1 − y n = δt 2 • Other Methods: Runge Kutta, Adams Bashforth, Backward differentiation, splitting The Heat Equation • Used to model diffusion of heat, species, • 1D • 2D • 3D ∂u ∂2u = ∂t ∂x 2 ∂u ∂2u ∂2u = + ∂t ∂x 2 ∂y 2 ∂u ∂2u ∂2u ∂2u = + + ∂t ∂x 2 ∂y 2 ∂z 2 • Not always a good model, since it has infinite speed of propagation • Strong coupling of all points in domain make it computationally intensive to solve in parallel ¨ Linear Schrodinger Equation • Used to model quantum mechanical phenomena and often appears in simplified wave propagation models • 1D i • 2D i • 3D i ∂u ∂2u = ∂t ∂x 2 ∂2u ∂2u ∂u = + ∂t ∂x 2 ∂y 2 ∂u ∂2u ∂2u ∂2u = + + ∂t ∂x 2 ∂y 2 ∂z 2 • Looks like a heat equation with imaginary time • Strong coupling of all points in domain make it computationally intensive to solve in parallel The Wave Equation • Used to model propagation of sound, light • 1D ∂2u ∂2u = ∂t 2 ∂x 2 • 2D • 3D ∂2u ∂2u ∂2u = + ∂t 2 ∂x 2 ∂y 2 ∂2u ∂2u ∂2u ∂2u = + + ∂t 2 ∂x 2 ∂y 2 ∂z 2 • Has finite speed of propagation • Finite signal propagation speed sometimes useful in parallelization since there is no coupling between grid points that are far apart, hence smaller communication requirements • If propagation speed is very fast, communication requirements still important The Poisson Equation • Used to model static deflection of a drumhead • 1D f (x) = • 2D f (x, y) = • 3D f (x, y, z) = ∂2u ∂x 2 ∂2u ∂2u + ∂x 2 ∂y 2 ∂2u ∂2u ∂2u + + ∂x 2 ∂y 2 ∂z 2 • No time dependence, but also arises in time discretizations of time dependent partial differential equations Burgers’ Equation • Simple model for gas dynamics, also traffic • Inviscid form (can develop shocks) ut = −uux • Viscous form ut = −uux + uxx Navier-Stokes Equations • Models motion of fluids • Consider incompressible case only (low speed motion relative to speed of sound) • ρ ∂u + u · ∇u = −∇p + µ∆u ∂t ∇ · u = 0. • p pressure, µ viscosity, ρ, density • 2D u(x, y ) = (u(x, y), v (x, y)) • 3D u(x, y , z) = (u(x, y, z), v (x, y, z), w(x, y, z)) Maxwells Equations • Model propagation of electromagnetic waves • ~t − ∇ × H ~ =0 D • ~t + ∇ × E ~ =0 B • with the relation between the electromagnetic components given by the constitutive relations: • ~ = εo (x, y, z)E ~ D • ~ = µo (x, y, z)H ~ B Summary on Models • Many diverse models • In 2D and 3D, parallel computing is very useful for getting numerical solutions in reasonable time • Typical applications are in physics, chemistry, biology and engineering • Related models also appear in social sciences, though usefulness for predictability of real world data is less clear Aim to simulate this, • https://www.flickr.com/photos/kitware/2293740417/in/pool-paraview/ • Visualization around Formula 1 Race Car by Renato N. Elias, Rio de Janeiro, Brazil or this • https://www.flickr.com/photos/kitware/2294528826/in/pool-paraview/ • Visualization around Formula 1 Race Car by Renato N. Elias, Rio de Janeiro, Brazil Maybe even this • https://youtu.be/FTlJ9l9KTkw • https://youtu.be/MsDgw8_cb90 • https://youtu.be/Su8qzC4HHVs A simple model: The Heat Equation The Heat Equation • How does the temperature of the rod vary with time? The Heat Equation: Model 1 • Middle of rod is initially hot due to previous heating (eg. oven in house that has just been turned off) • Temperature on the two sides is 0 (winter and cold outside the house) • Assume discrete uniformly space time, and discrete space with molecules at each coordinate point • Assume temperature in bar is directly proportional to heat, with the same proportionality constant throughout the bar • At next time step, a molecule will • with probability 0.5 not transfer its current heat • with probability 0.25 transfer half of its heat to its left neighbor • with probability 0.25 transfer half of its heat to its right neighbor The Heat Equation: Model 1 • Let Tin denote the temperature at position i at time n. Particles at the end have a fixed temperature/heat they can transfer, but they always remain at the same temperature. n+1 Ti = Not transfer heat 0.0Ti−1 + Ti + 0.5Ti+1 n n probability 0.75 × 0.5 × 0.25 n 0.5Ti−1 n n probability 0.25 × 0.5 × 0.75 + n Ti + 0.0Ti+1 n n n probability 0.25 × 0.5 × 0.25 n n n probability 0.75 × 0.5 × 0.75 0.5Ti−1 + Ti + 0.5Ti+1 0.0Ti−1 + Ti + 0.0Ti+1 Transfer heat left n n n probability 0.75 × 0.25 × 0.25 n n n probability 0.25 × 0.25 × 0.75 0.0Ti−1 + 0.5Ti + 0.5Ti+1 0.5Ti−1 + 0.5Ti + 0.0Ti+1 n 0.5Ti probability 0.25 × 0.25 × 0.25 n 0.5Ti−1 + + n 0.5Ti+1 n n probability 0.75 × 0.25 × 0.75 0.0Ti−1 + 0.5Ti + 0.5Ti+1 n n probability 0.75 × 0.25 × 0.25 n 0.5Ti−1 n probability 0.25 × 0.25 × 0.75 n 0.0Ti−1 + 0.5Ti + 0.0Ti+1 Transfer heat right n + n 0.5Ti + 0.0Ti+1 n n n probability 0.25 × 0.25 × 0.25 n n n probability 0.75 × 0.25 × 0.75 0.5Ti−1 + 0.5Ti + 0.5Ti+1 0.0Ti−1 + 0.5Ti + 0.0Ti+1 The Heat Equation: Model 1 • DEMO The Heat Equation: Model 2 • Probability is hard! • It takes a long time • May only interested in average temperature/heat, not instantaneous fluctuations • Take averages n Tin+1 = (0.25 × 0.5)Ti−1 + (0.5 + 0.25 × 0.5 + 0.25 × 0.5)Tin n + (0.25 × 0.5)Ti+1 n n Tin+1 = 0.125Ti−1 + 0.75Tin + 0.125Ti+1 or Tin+1 = Tin + n − 2T n + T n Ti−1 i i+1 8 • The last line looks like a difference approximation for a differential equation The Heat Equation: Model 3 • Let us find a differential equation! • Make the time increment small n − 2T n + T n Ti−1 i i+1 8 n − 2T n + T n Ti−1 Tin+1 − Tin i i+1 = δt 8δt Tin+1 = Tin + The Heat Equation: Model 3 • Let us find a differential equation! • Make the space increment small n T n − 2Tin + Ti+1 Tin+1 − Tin = i−1 δt 8δt Tin+1 − Tin (δx)2 = δt 8δt n −T n Ti−1 i δx − δx n Tin −Ti+1 δx • Let δx → 0 and δt → 0 such that (δx)2 = δt we get ∂T 1 ∂2T = ∂t 8 ∂x 2 • Constant in continuum formulation depends on physics and is usually measured experimentally, or determined from microscopic (many particle) simulations. The Heat Equation: Numerical Solution Methods • Probabalistic • Finite Difference • Finite Volume • Finite Element • Spectral Finite Difference Method • Approximate derivatives by difference quotients • Simple method to derive and implement • Convergence rates tend not to be great • Difficult to use for complicated geometries • Tends to scale well since communication requirements are low Finite Difference Method for Heat Equation • ut = 8uxx • Using backward Euler time stepping: n+1 u n+1 − 2uin+1 + ui+1 uin+1 − uin = 8 i−1 δt (δx)2 • Using forward Euler time stepping (strong stability restrictions): n u n − 2uin + ui+1 uin+1 − uin = 8 i−1 δt (δx)2 Finite Difference Method for Heat Equation • Simple method to derive and implement • Hardest part for implicit schemes is solution of resulting linear system of equations • Explicit schemes typically have stability restrictions or can always be unstable • Convergence rates tend not to be great – to get an accurate solution, a large number of grid points are needed • Difficult to use for complicated geometries • Tends to scale well since communication requirements are low Finite Volume Method for Heat Equation 2 ∂u ∂v • = ∂∂xu2 or ∂u ∂x = v and ∂t = ∂x • Consider a cell averaged integral, then use implicit ∂u ∂t midpoint rule R xi+1 Rx ux dx = xi i+1 v dx Rx ui+1 − ui = xi i+1 v dx xi n+1 n −u n+1 −u n ui+1 +ui+1 i i 2 = n+1 n +v n+1 +v n vi+1 +vi+1 i i δx 2 Rx ut dx = xi i+1 vx dx xi R xi+1 ut dx = vi+1 − vi xi R xi+1 n+1 n −u n ui+1 +uin+1 −ui+1 i δx 2δt = n+1 n −v n+1 −v n vi+1 +vi+1 i i 2 • Several ways of approximating the integrals. The one above is a little unusual, most finite volume schemes use left sided or right sided approximations. Finite Volume Method for Heat Equation • For implicit schemes, hardest part is solving the system of equations that results • Explicit schemes parallelize very well, however a large number of grid points are usually needed to get accurate results • Automated construction of simple finite volume schemes is possible, making them popular in packages • No convergence theory for high order finite volume schemes • Tricky to do complicated geometries accurately Finite Element Method for Heat Equation ∂2u • ∂u ∂t = ∂x 2 P ˜ (x, t) = i φi (x)Ti (t), where φi (x) is • Assume u(x, t) ≈ u zero for x > xi+1 and x < xi−1 . • A simple choice is ψ is the triangle hat function, xi+1 − x for x ∈ [xi , xi+1 ] and x − xi−1 for x ∈ [xi−1 , xi ] • Now need to find Ti (t), which will be evaluated using finite differences. Finite Element Method for Heat Equation • Consider multiplying the heat equation by a polynomial, φj that is only non zero over a few grid points, then integrate by parts, R xi+1 ∂ 2 u˜ ˜ ∂u xi−1 ∂t φi dx = xi−1 ∂x 2 φi dx R xi+1 ∂ u˜ R xi+1 ∂ u˜ ∂φi xi−1 ∂t φi dx = − xi−1 ∂x ∂x dx R xi+1 Finite Element Method for Heat Equation • Then use implicit midpoint rule R xi+1 R xi+1 ∂ u˜ ∂φi = − xi−1 ∂x ∂x dx R xi+1 1 ∂ u˜n ∂ u˜n+1 ∂φi ˜ n+1 −u ˜n u φ dx = − i δt ∂x + ∂x ∂x dx xi−1 2 ˜ ∂u xi−1 ∂t φi dx R xi+1 xi−1 1 2 n+1 ˜ n ˜i+1 u −ui+1 φi |xi+1 δt n+1 ˜ n ˜i−1 u −ui−1 φi |xi−1 δt 2δx ∂φi ˜ n+1 ˜n ∂u 1 1 ∂u + ∂x 2δx = −2 2 ∂x x ∂x x xi+1 i+1 i+1 n+1 ∂φi ˜n u − 21 12 ∂ u˜∂x + ∂∂x 2δx ∂x xi−1 + xi−1 xi−1 • Since φ are known before hand, can re-write this as matrix vector products, so need to solve a linear system at each time step. Finite Element Method for Heat Equation • Several other ways of approximating the integrals, can extend to multiple dimensions. • Weak formulation allows for solution of equations where second derivative is not naturally defined • Large mathematical community developing convergence theory for these methods • Well suited to complicated geometries • Rather difficult to implement compared to other schemes because of integrals that need to be computed • Used in many codes, but typically codes are hand written to obtain high efficiency • For implicit time discretizations, solving the linear system of equations that results can be most time consuming part Fourier Spectral Method for Heat Equation 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 • ∂2u ∂u = ∂t ∂x 2 • 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 exp(ikx). = (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 Key Concepts • Differential equations used to model physical phenomena - approximate model • Microscopic model can suggest ways to find approximate solutions on computers • Several methods introduced, probabilistic, finite difference, finite element, spectral (Fourier transform) • Need to check simulations for discretization errors, then before using for prediction, understand modeling errors • Time stepping, Forward Euler, Backward Euler, Implicit Midpoint Rule, Implicit-Explicit References • Barth, T. and Ohlberger, M. “Finite Volume Methods: Foundation and Analysis” https://archive.org/details/nasa_techdoc_20030020790 • Boyd, J.P. “Chebyshev and Fourier Spectral Methods” 2nd. edition http: //www-personal.umich.edu/˜jpboyd/BOOK_Spectral2000.html • 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/ Cooley J.W. and Tukey J.W. “An algorithm for the machine calculation of complex Fourier series” Math. Comput. 19, 297–301, (1965) ¨ T. and Herbin, R. “Finite Volume Methods” Eymard, R. Gallouet, http://www.cmi.univ-mrs.fr/˜herbin/BOOK/bookevol.pdf Greenbaum A. and Chartier T.P. Numerical Methods Princeton University Press (2012) 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) Ketcheson, D. “Hyperpython” https://github.com/ketch/HyperPython/ Lawler, G.F. “Random Walk and the Heat Equation” American Mathematical Society (2010) Macqueron, C. “Computational Fluid Dynamics Modeling of a wood-burning stove-heated sauna using NIST’s Fire Dynamics Simulator” http://arxiv.org/abs/1404.6774 References • McGrattan, K. et al. Fire Dynamics Simulator http://www.nist.gov/el/fire_research/fds_smokeview.cfm • Petersen W.P. and Arbenz P. Introduction to Parallel Computing Oxford University Press (2004) • Quarteroni, A. Sacco, R. and Saleri, F. “Numerical Mathematics” http://link.springer.com/book/10.1007%2Fb98885 • Sauer T. Numerical Analysis Pearson (2012) • Suli, ¨ E. “Finite element methods for partial differential equations” people.maths.ox.ac.uk/suli/fem.pdf • Trefethen, L.N. “Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations” http://people.maths.ox.ac.uk/trefethen/pdetext.html • Vainikko E. “Scientific Computing Lectures” https://courses.cs.ut.ee/2013/scicomp/fall/Main/Lectures ¨ • Vainikko E. Fortran 95 Ja MPI Tartu Ulikool Kirjastus (2004) http://kodu.ut.ee/˜eero/PC/F95jaMPI.pdf Acknowledgements • Oleg Batra˘sev • 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 Extra Slides 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) 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 =65536 # 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 / numpoints I2dsquare = I2dsquare / numpoints E s t i m E r r o r =4.0∗numpy . s q r t ( ( I2dsquare − I 2 d ∗∗2)/ numpoints ) # e s t i m a t e d e r r o r I 2 d = I 2 d ∗4.0 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 1: 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 37.09 44.49 44.50 43.80 44 Error Estimate +/- 10.89 +/- 2.40 +/- 0.60 +/- 0.14 0.0 See http://en.wikibooks.org/wiki/Parallel_ Spectral_Numerical_Methods/Introduction_to_ Parallel_Programming for more information
© Copyright 2024