Description

Lecture 10
PDE and Programming – 1
Jung-Il Choi

Lecture 10
• Partial Differential Equation • 1D PDEs
 1D Heat equation
 Semi-discretization
 Stability analysis
 Eigenvalue analysis
 Modified wavenumber analysis
 von Neumann analysis
 Accuracy via modified equation
 Example 1
 1D Wave equation
 Semi-discretization
 Stability analysis
 Modified wavenumber analysis
 Example 2
Contents
Lecture 11
• Multi-dimension
 Heat equation
 Implicit methods in higher
 Approximate factorization
 Stability analysis
 Alternating direction implicit methods (ADI)
 Poisson equation
 Iterative solution methods
 Point Jacobi method
 Gauss-Seidel method
 Successive over relaxation method (SOR)
 Non-linear PDEs
2 / 41
• Partial Differential Equation (PDE)
An equation stating a relationship between a function of two or more independent variables and the partial derivatives of this function with respect to these independent variables.Non linear의 대표적인 예
0: 2D Laplace equation
𝑢 0
 𝑓 𝑥, 𝑦: 2D Poisson equation
 → 𝑢 𝑢 𝑓
 𝑐 : 1D Diffusion equation
 (→ 𝑢 𝑐 𝑢 )
 𝑐 : 1D Wave equation
 (→ 𝑢 𝑐 𝑢 ) 이 경우는 일정한 속도와 방향으로 움직임 속도가 빠르면 빠르게, Diffusion equation 느리면 느리게
하는 wave equation
ΔΦ 4𝜋𝐺𝜌
https://en.wilipedia.org/wiki/
3 / 41

• Solution and linearity(or nonlinearity) of PDEs
 The solution of a PDE in some region 𝑅 of the domain of interest, 𝐷 𝑥⃗, 𝑡
 the particular function, 𝑢 𝑥⃗, 𝑡 satisfies the PDE in 𝑅,
 and the initial and/or boundary conditions specified on the boundaries of 𝑅 ⊂ 𝐷.
평균에 대해서는 보장을 받을 수 있음 (analysis가능)
initial을 무조건 정.확.히 알아야 한다는 것이 아님! 여러번의 trial을 해야함
 Linear PDE https://m.blog.naver.com/pmw9440/221442252220
 All partial derivatives appear in a linear form (first degree in the unknown function 𝑢 and its derivatives)
 “AND” none of the coefficients depend on the dependent variable
𝑢 𝑢 0, 𝑢 c 𝑢 , 𝑎𝑢 𝑏𝑥𝑢 0
 Nonlinear PDE
Variable coefficient linear PDE
 The derivatives appear in a nonlinear form depends on “independent” variable
 “OR” the coefficients depend on the dependent variable
𝑢𝑢 𝑏𝑢 0, 𝑎𝑢 𝑏𝑢
4 / 41
• Order of PDE
The highest-order derivative
𝜕 𝑢 𝜕𝑢
0 → 2 𝑜𝑟𝑑𝑒𝑟 𝜕𝑥 𝜕𝑦
• Homogeneous vs Nonhomogeneous
Homogeneous PDE
each of the terms contains 𝑢 or the dependent variables or its partial derivatives.
𝑢 𝑢 0
𝑢𝑢 𝑏𝑢 0 𝑎𝑢 𝑏𝑢 0
Nonhomogeneous PDE
∇ 𝑢 𝑢 𝑢 𝑢 𝑓 𝑥, 𝑦, 𝑧
5 / 41
• Classification of PDEs using characteristics analysis
Consider the general quasilinear 2nd order nonhomogeneous PDE in 2D:
𝐴𝑢 𝐵𝑢 𝐶𝑢 𝐹 𝑥, 𝑦, 𝑢, 𝑢 , 𝑢
 Parabolic equation (𝐵
 Hyperbolic equation (𝐵 4𝐴𝐶

4𝐴𝐶 Quasilinear : linear in the highest-order derivative
0) 에너지가 많다면, 온 사방에 균일하게 나눠줌모래성이 사라지는 그림 생각하면 좋아 모양이 포물선 같아서!
𝑐 ∇ 𝑢 : Diffusion equation 0)
파도가 쳐서 오는 모양!
u(x,t)을 y라는 식을 써서 평행이동 해버리자
𝑐 ∇ 𝑢 : Wave equation t축과 x축이 있을 때, 매개체가 c:속도
c에 t를 곱하면 간 거리가 되니, y = x – ct로 하면 1변수로 바뀜
u(x,t) ==> f(x-ct) + f(x+ct) : solution이 dependent함
Elliptic equations (𝐵 4𝐴𝐶 0)
∇ 𝑢 0
∇ 𝑢 𝑓 𝑥⃗ : Laplace equation (homogeneous) and Poisson equation (nonhomogeneous)
Heat equation의 steady stationary condition이면 laplace eq.
—> 시간에 대한 요소가 없음
—> diffusion eq. solution이 변하지 않을 때까지 가는 것 ( time scale이 없어
6 / 41
• Classification of PDEs using characteristics analysis
The terminology elliptic, parabolic, and hyperbolic chosen to classify PDEs reflects the analogy between the from of the discriminant 𝑩𝟐 𝟒𝑨𝑪
(from the idea of d’Alembert’s solution, methods of characteristics) And that which classifies conic section.
𝐴𝑥 𝐵𝑥𝑦 𝐶𝑦 𝐷𝑥 𝐸𝑦 𝐹 0
Type Defining condition Examples
Parabolic 𝐵 4𝐴𝐶 0 Diffusion equation
Hyperbolic 𝐵 4𝐴𝐶 0 Wave equation
Elliptic 𝐵 4𝐴𝐶 0 Laplace/Poisson equation
𝐴𝑢 𝐵𝑢 𝐶𝑢 𝐹 𝑥, 𝑦, 𝑢, 𝑢 , 𝑢
https://en.wilipedia.org/wiki/
7 / 41
• Classification of PDEs using characteristics analysis
 Characteristics
 Propagate behavior of each fixed point on the space at the “Hyper” space( 𝑛 1 D space for 𝑛D PDE) Information, 𝑢 (velocity, temperature, pressure etc.) propagates along path.
 Are there any points in the solution domain 𝐷 𝑥, 𝑦 passing through a general point 𝑃 along which the second derivatives of 𝑢 x, y are multivalued or discontinuous (kernel space)?
 Homogeneous solution
 If there are such paths, they are called path of information propagation Or Characteristics
8 / 41
• Classification of PDEs using characteristics analysis
Chain rule & Homogeneous solution (Kernel space)
𝑑 𝑢 𝑢 𝑑𝑥 𝑢 𝑑𝑦
𝐴
𝑑𝑥
𝐵
𝑑𝑦
𝑑𝑥 𝐶

𝑑𝑦 𝑢
𝑢
𝑢
𝑢 𝑑𝑦
𝐴 𝐵 𝐶
⇒ det 𝑑𝑥 𝑑𝑦 0
𝑑𝑥 𝑑𝑦
0
Discriminant Characteristics Type
𝐵 4𝐴𝐶 0 Real & Repeated Parabolic
𝐵 4𝐴𝐶 0 Real & Distinct Hyperbolic
𝐵 4𝐴𝐶 0 Complex Elliptic
𝑑𝑦 𝐵 𝐵 4𝐴𝐶
⇒
𝑑𝑥 2𝐴
9 / 41
• Classification of PDEs using characteristics analysis
 Parabolic PDEs have one real repeated characteristic path (Critical damping, diffusing)
 Hyperbolic PDEs of two real distinct characteristic paths (Overdamping, diffusing)
 Elliptic PDEs have no real characteristic paths (Oscillatory)
10 / 41

One-dimensional PDEs

• 1D Heat equation
 Semi-discretization
 Temporal discretization
 Stability analysis
 Eigenvalue/Eigenvector analysis
 Modified wavenumber analysis
 von Neumann analysis
 Accuracy via modified equation
 Example 1
• 1D Wave equation
 Semi-discretization
 Stability analysis
Modified wavenumber analysis
 Example 2
11 / 41
Full로 하려니 너무 써야하는게 많아..!! 1D Heat equation
• Semi-discretization – Solving a PDE as a system of ODEs
 Numerical methods for PDEs are straightforward extensions of methods developed for initial and boundary value problems in ODEs.
 That is, a PDE can be converted to a system of ODEs by using finite difference methods for the derivatives in all but one of dimensions.
 Consider the one-dimensional diffusion(or heat equation)
𝜕𝜙 𝜕 𝜙 Initial condition : 𝜙 𝑥, 0 𝑔 𝑥

𝜕𝑡 𝛼 𝜕𝑥 Boundary condition : 𝜙 0, 𝑡 𝜙 𝐿, 𝑡 0
 Discretization of the Domain with 𝑁 intervals → 𝑁 1 uniformly spaced grid points
Δ𝑥 Δ𝑥

𝑥 𝑥 𝑥 𝑥 𝑥 𝑥 𝑥 𝑥 𝑥 Δ𝑥
𝑗 0, 𝑗 𝑁 are the boundaries
𝑗 1, 2, 3, ⋯ , 𝑁 1 are interior points
12 / 41

공간차분 먼저 (x(j)시점)
Semi-discretization
 Let’s use the second-order central difference scheme to the second derivative.
𝑑𝜙 𝜙 2𝜙 𝜙
𝛼 , 𝑗 1,2,3, ⋯ , 𝑁 1 𝑑𝑡 Δ𝑥
Where 𝜙 𝜙 𝑥, 𝑡
 A system of 𝑁 1 ordinary differential equations
 Space derivatives for fixed time (→ Semi-discretization) and solving time marching as solving ODEs.
 Can be written in matrix form as:
𝑑𝜙
𝐴𝜙
𝑑𝑡
 Where as 𝜙 are the (time-dependent) elements of the vector 𝜙, and 𝐴 is an 𝑁 1 𝑁 1 tridiagonal
matrix.
13 / 41
Semi-discretization
2 1 ⋯
𝛼 1 2 1
𝐴
Δ𝑥 ⋱ ⋱ ⋱
1 2
→ 𝑁 1 𝑁 1 tridiagonal matrix which is symmetric
 The result is a system of ODEs that can be solved using any of the numerical methods introduced for ODEs, such as Euler methods, RK formulas or multi-step methods.
 However, when dealing with systems, we should be concerned about stability.
 The range of the eigenvalues of 𝐴 determines whether the system is stable.
𝐴𝜙 𝜆𝜙 → det 𝐴 𝜆𝐼 0
14 / 41
Temporal discretization
• (Recall) Various time advancement schemes
𝑑𝜙
𝐴𝜙
𝑑𝑡
 Forward Euler scheme
𝜙 𝜙
𝐴𝜙
Δ𝑡
 Backward Euler scheme
𝜙 𝜙
𝐴𝜙
Δ𝑡
 Crank-Nicolson scheme
𝜙 𝜙
𝐴
Δ𝑡
15 / 41

• Eigenvalue/Eigenvector analysis
 (Recall) Diagonalization, Eigenvalues, Λ & Eigenvectors(Eigenfunctions), 𝑋
Diagonalization (Decoupling)
Suppose 𝐴 has the eigenvalues (𝑖𝑒. 𝐴 is diagonalizable),
𝑋 𝐴𝑋 Λ → 𝐴 𝑋Λ𝑋
𝑑𝜙 𝑑𝜙𝑑𝜙
𝐴𝜙 ⇒ 𝑋Λ𝑋 𝜙 ⇒ 𝑋 Λ𝑋 𝜙
𝑑𝑡 𝑑𝑡𝑑𝑡
𝑑 𝑋 𝜙𝑑𝜓
Λ𝑋 𝜙 ⇒ Λ𝜓, 𝜓 𝑋 𝜙
𝑑𝑡𝑑𝑡
𝑑𝜓
𝜆𝜓 ⇒ 𝜓 𝑐𝑒 ⇒ 𝜓 𝑐
𝑑𝑡
𝑒
0
⋮ 𝑐 0
𝑒
⋮ ⋯ 𝑐 0
0
⋮
0 0 𝑒
𝜓 𝑐 𝑒 → 𝜙 𝑋𝜓
16 / 41
 Analytical expressions of eigenvalues of the matrix 𝐴
𝜋𝑗
𝜆2 2 cos, 𝑗 1,2,3, ⋯ , 𝑁 1
𝑁
 The eigenvalue with the smallest (𝑗 1) and the largest magnitude (𝑗 𝑁 1) is:
𝛼 𝜋
𝜆2 2 cos, 𝜆2 2 cos Δ𝑥 𝑁

 For large, 𝑁, the Taylor series expansion for cos converges rapidly, and cos converges to -1.
𝜋 1𝜋 𝑁 1
cos 1 ⋯ , cos 𝑐𝑜𝑠𝜋 1
𝑁 2!𝑁
 Using the first two terms in the expansion then,
𝜋 𝛼
𝜆 ,
𝑁 Δ𝑥
17 / 41
𝜋 𝛼
𝜆 4 ,
𝑁 Δ𝑥 𝜆 𝛼
4
Δ𝑥
 The ratio of the eigenvalue with the largest modulus to that with the smallest modulus is :

 For large N, the system is unstable!
18 / 41
frequency : 1초당 진동수 Wavenumber : wave의 숫자
• Modified wavenumber analysis 공간차분을 어떻게 하냐에 대해서 결정됨
Let revisit the heat equation,
𝛼
cos 𝑘Δ𝑥
𝑑𝑡
𝑑𝜓
𝛼𝑘′ 𝜓
𝑑𝑡
2
𝑘′1 cos 𝑘Δ𝑥
Δ𝑥
 𝛼𝑘 𝜆
𝜓 𝜆𝜓
 Using the forward Euler for time advancement,
2𝛼
1 cos 𝑘Δ𝑥 𝜓
Δ𝑡 Δ𝑥
2 2
Δ𝑡 ⇒ Δ𝑡
𝜆 2𝛼 1 cos 𝑘Δ𝑥
Δ𝑥
 Since, 1 cos 𝑘Δ𝑥 1, the worst-case scenario is :
Δ𝑥
Δ𝑡
2𝛼
20 / 41
𝛼𝑘
Δ𝑡
1 𝜓
 For the stability analysis,
𝑑𝜓
𝛼𝑘′ 𝜓
𝑑𝑡
2
𝑘′1 cos 𝑘Δ𝑥
Δ𝑥
𝜓 𝜎𝜓
Where 𝜎 ⇒ 𝜎 1
Crank-Nicolson is unconditionally stable
21 / 41
𝑑𝜓
𝛼𝑘′ 𝜓
𝑑𝑡
2
𝑘′1 cos 𝑘Δ𝑥
Δ𝑥
 Using backward Euler
𝛼𝑘 𝜓
Δ𝑡
1 𝜓 𝜓
 For the stability analysis,
𝜓 𝛾𝜓
1 1
Where 𝛾 ⇒ 𝛾 1
Backward Euler is unconditionally stable
 However, in contrast to Crank-Nicolson, 𝜎 → 0 when Δ𝑡 → ∞. That is, the solution does not exhibit undesirable oscillations (although it would be inaccurate).
22 / 41
Consider the wave equation,
𝜕𝑢 𝜕𝑢
𝑐 0, 0 𝑥 𝐿, 𝑡 0
𝜕𝑡 𝜕𝑥
 Assuming, 𝑢 𝑥, 𝑡 𝑣 𝑡 𝑒
𝑑𝑣𝑑𝑣
𝑒 𝑖𝑘𝑐 𝑒 𝑣 ⇒ 𝑖𝑘𝑐 𝑣 Modified wavenumber
𝑑𝑡𝑑𝑡
 Semi-discretized equation with central difference scheme,
𝑑𝑢 𝑢 𝑢 𝑑𝑣 sin 𝑘Δ𝑥
𝑐 0 ⇒ 𝑖𝑐 𝑣 𝑖𝑐𝑘𝑣
𝑑𝑡 2Δ𝑥 𝑑𝑡 Δ𝑥
23 / 41
• von Neumann stability analysis
 Matrix stability analysis using the eigenvalues of the matrix obtained from a semi-discretization of
PDE
 This is only available for very simple matrices Consider full discretization of PDE
 von Neumann stability analysis does not account for the effect of boundary conditions; periodic boundary conditions are assumed.
 Linear, constant coefficient differential equations with uniformly spaced spatial grids.
𝜕𝜙 𝜕 𝜙
𝛼
𝜕𝑡 𝜕𝑥
 Second-order central difference with the explicit Euler method
𝜙 𝜙 𝜙 2𝜙 𝜙
𝛼
Δ𝑡 Δ𝑥
 Assuming a solution of the form
𝜙 𝜎 𝑒
24 / 41

𝜙𝜙2𝜙𝜙
𝜎 𝑒 𝜎 𝑒 2𝑒 𝑒
Where 𝑥 𝑥 Δ𝑥and 𝑥 𝑥 Δ𝑥.
𝛼Δ𝑡
𝜎 12 cos 𝑘Δ𝑥 2
Δ𝑥
For stability, 𝜎 1
2 cos 𝑘Δ𝑥 1
𝛼Δ𝑡Δ𝑥
2 cos 𝑘Δ𝑥 2 → Δ𝑡
Δ𝑥
The worst (or the most restrictive) case occurs when cos 𝑘Δ𝑥
Δ𝑥
Δ𝑡
2𝛼
25 / 41

Accuracy via modified equation
 By 1 𝛼 2 ,
𝜕 𝜙
𝛼
𝜕𝑥
 Consider the (unsteady) heat equation (or 1D diffusion equation) given by:
 𝑥 𝑥 Δ𝑥
 Discretized equation:

28 / 41

 Discretized equation:
 The PDE has been converted to a system of ODEs
 Time advancement
 Using forward Euler,
𝑇 𝑇 Δ𝑡𝐹 𝑇 , 𝑡
1,2,3, ⋯ , 𝑁 1
𝛼
𝑇 2𝑇
Δ𝑥
𝛼
𝑇 2𝑇 𝑇
Δ𝑥
⋮
𝜋 1 𝑒 sin 𝜋𝑥
𝜋 1 𝑒 sin 𝜋𝑥
29 / 41

The stability of the numerical solution for time advancement depends on the eigenvalue of the system having the largest magnitude:
𝛼
𝜆 4
Δ𝑥
 When forward Euler is used for real and negative 𝜆:
2 Δ𝑥
Δ𝑡
2𝛼
 For 𝛼 1 and Δ𝑥 0.05,
Δ𝑡 0.00125
30 / 41
• Pseudo-code & Code

31 / 41
• result
 Δ𝑡 0.001, Δ𝑥 0.05, 𝛼 1

 Δ𝑡 0.0015, Δ𝑥 0.05, 𝛼 1

32 / 41

(dissipatio i 에너지가점점없어집 n dispersion: 에너지가점점많아짐? p 시간이지나도모양이변하지않아야함 .
 Consider a semi-discretization of the following first-order wav1diD Wave equationssipation되는卜挑昨 e equation :차분화할때 ,없내그을래가수원프도하모양이있는격잖아자위?!치에
 aka the convection/transport equation 그래서어렵다. (FDM),
𝜕𝑢 𝜕𝑢 Initial condition : 𝑢 𝑥, 0 𝑓 𝑥
𝑐 0
𝜕𝑡 𝜕𝑥 Boundary condition : 𝑢 0, 𝑡 𝜙 𝐿, 𝑡 0 岬州싶은데’
 A simple model equation for the convection phenomena. l l h t l l l r
 The exact solution is such that an initial disturbance in the domain (𝑢 𝑥, 0 ) simply propagates with the constant convection speed 𝑐 in the positive (or negative) 𝑥-direction.
33 / 41
𝜕𝑢 𝜕𝑢
𝑐 0
𝜕𝑡 𝜕𝑥
Semi-discretization
Semi-discretization
 Assume 𝑐 0, and using central difference scheme, 시간에대한차분X c 공간에대한차분0.
𝑑𝑢 𝑢 𝑢
𝑐 0
𝑑𝑡 2Δ𝑥 의 System o f ODE가됨 .
 In matrix from,
𝑑𝑢
𝐴𝑢
𝑑𝑡
0 1 ⋯
1 0 1 where 𝐴 𝑁 1 𝑁 1 tridiagonal matrix which is not symmetric
⋱ ⋱ ⋱
1 0
 From analytical consideration, no boundary condition is prescribed at 𝑥 𝐿.
 However, a special numerical boundary treatment is required at 𝑥 𝐿 owing to the use of central differencing in this problem.
 A typical well-behaved numerical boundary treatment at 𝑥 𝐿 slightly modifies the last row of the coefficient matrix 𝐴, but we will ignore it for now.
34 / 41

 Thus, the eigenvalues of the matrix resulting from semi-discretization of the convection equation are purely imaginary N 바람이부는방향으로N . balanced된 stability써도
upwind 좋지않아서결국 Upwind .
𝑐 𝜋𝑗
𝜆 𝑖𝜔, where 𝜔 cos
Δ𝑥 𝑁
 The solution is a superposition of modes, where each mode’s temporal behavior is given by 𝑒 Oscillatory or sinusoidal(non-decaying) character.
35 / 41
𝑑𝑣
𝑖𝑐𝑘𝑣
𝑑𝑡
sin 𝑘Δ𝑥
𝑘
Δ𝑥
Stability analysis
 Leap flog method for time advancement,
𝑣 𝑣 sin 𝑘Δ𝑥

 Consider the numerical solution to the homogeneous convection equation
𝜕𝑢 𝜕𝑢
𝑐 0
𝜕𝑡 𝜕𝑥 0
𝑡 𝑥 𝐿 0
 Initial conditions: 𝑢 𝑥, 0 𝑒 .
 Boundary conditions: 𝑢 0, 𝑡 0
 Although the proper spatial domain for this PDE is semi-infinite, numerical implementation requires a finite domain.
 Thus, we arbitrarily truncate the domain to 0 𝑥 1
 Semi-discretized equation using a 2nd order central difference scheme:
𝑑𝑢 𝑢 𝑢
𝑐 0
𝑑𝑡 2Δ𝑥
37 / 41
• Pseudo-code & Code

𝑥
𝑡
𝑢 u ← e
0
𝑐 ← 1
𝑑𝑢𝑑𝑡
𝑑𝑢𝑑𝑡
𝑑𝑢𝑑𝑡 Program solve_wave_eq
← 0, 𝑥 ← 0.75
← 0, 𝑛𝑡 ← 21 𝑑𝑡 ← 0.01, 𝑑𝑥 ← 0.01
.
← 0
𝐶𝑎𝑙𝑙 𝐸𝑢𝑙𝑒𝑟 𝑜𝑟 𝑅𝐾4 𝑡, 𝑢, 𝑛𝑡, 𝑑𝑡, 𝑤𝑎𝑣𝑒𝑒𝑞, 𝑥, 𝑑𝑥
End program
Function waveeq(T, t, x, dx)
←
←
𝑟𝑒𝑡𝑢𝑟𝑛 𝑑𝑢𝑑𝑡
End function

38 / 41

𝑥
𝑡
𝑢 u ← e
0
𝑐 ← 1
𝑑𝑢𝑑𝑡
𝑑𝑢𝑑𝑡
𝑑𝑢𝑑𝑡 Program solve_wave_eq
← 0, 𝑥 ← 0.75
← 0, 𝑛𝑡 ← 21 𝑑𝑡 ← 0.01, 𝑑𝑥 ← 0.01
.
← 0
𝐶𝑎𝑙𝑙 𝐸𝑢𝑙𝑒𝑟 𝑜𝑟 𝑅𝐾4 𝑡, 𝑢, 𝑛𝑡, 𝑑𝑡, 𝑤𝑎𝑣𝑒𝑒𝑞, 𝑥, 𝑑𝑥
End program
Function waveeq(T, t, x, dx)
←
←
𝑟𝑒𝑡𝑢𝑟𝑛 𝑑𝑢𝑑𝑡
End function

• Pseudo-code & Code
39 / 41
𝑓𝑜𝑟𝑤𝑎𝑟𝑑 𝐸𝑢𝑙𝑒𝑟: 𝐶𝐹𝐿 1
𝑅𝐾4: 𝐶𝐹𝐿 2.83
• result Δ𝑡 0.01, Δ𝑥 0.01, 𝑐 1)
Euler method Fourth order Runge-Kutta method

Q&A Thanks for listening