Hyperbolic Problems

Numerical methods for differential equations

Hyperbolic PDEs

Hyperbolic partial differential equations (PDEs) have a totally different nature with respect to others. Let's start our analysis from conservation laws:

_∂u/_∂t + _∂(f(u))/_∂x = 0

^∂/^∂t > c₁ > c₂ > c_b f(x)

_∂g(u)/_∂x generally could be q f(u,x)

_∂u/_∂t + a(u) _∂u/_∂x conservative form

If a(u) = const., then:

• du/dt + a du/dx = 0
• f(u) = 1/2 u² ⇒ Burgers
• _∂u/_∂t + 2u _∂u/_∂x

a(u) = x₀ ⇒ dx/dt

x(t) - x₀ = a t ⇒ dx/dt = a t ⇒ Sub = x₀

Periodic boundary conditions: I trace back to the foot

Generic (non-periodic) boundary conditions: only on inflow

➔ Inflow is where a·n

Exist Γ such that A:TΛT^-1

Condition for system to be hyperbolic if A can be reduced to diag through Γ

Non-conservative form

∂u/∂t + j(u)/∂x = 0

Generally could be f(u,x)

∂u/∂t + ∂(a(u)u)/∂x = 0 ⇒ ∂u/∂t + a(u) ∂u/∂x

If a(u) = const., then:

• ∂u/∂t + a ∂u/∂x = 0
• If f(u) = 1/2 u² ⇒ Burgers

∂u/∂t + ∂u/∂t + a(u) ∂u/∂x

∂(u)/∂t; a(u) = u ⇒ ∂u/∂t + a(u) ∂u/∂x = 0

∂u/∂t + a(u) ∂u/∂x = 0 ⇒ total derivative is zero!

x(t) - x₀ = a(u) t ⇒ x(t) = x₀ + a(u) t

Characteristic lines: x(x) = x₀

Unbounded, i.e. u(x,0) = u₀(x)

∂u/∂t + a ∂u/∂x = 0

u(x,t) = u₀(x-at)

Periodic b.c.: I trace back to the foot

Generic (non-periodic) b.c.: only on inflow

Foot of characteristic line from (x,t) ⇒ (x-at,0)

INFLOW is where aⁿ < 0

∂Ω

Going to multi-dimensional case, ∂u/∂t + A ∂u/∂x = 0, A(p×p) u,p-vector

Exist ∇ such that A=TΛT^-1 where Λ is a diagonal matrix containing all eigenvalues of A, eigenvalues being real numbers Λ = diag(λ₁, ..., λ_p)

Condition for system to be hyperbolic is A can be reduced to diag through T.^T_i^-1

T^-1 ∂^{u/ _∂t} + ΛΓ^-1 ∂^{u/ _∂x} = 0

w = Γ^-1u define ^{w/ _∂t} + Λ ^{w/_{∂ x}} = 0

-> It's a decoupled system if Λ-∇ΓΛΓ^-1 holds

New variables w are called characteristic variables

1D characteristic lines Suppose p=3, p p system

• λ >0
• λ >0
• λ <0

Boundary conditions must be provided on proper points (inflow & outflow)

D(x,t)={x-λ_it, i=1,...,p}

w = Γ^-1 u ; A = ∇ΓΛ^-1 -> λ_i ∂^{w_i/ _∂x} = ∂^{w_i/ _∂t}

If A is diag, Γ will be a linear combination, Γ and T^-1 will be just coefficients, related to eigenvalues of A

If A(u) in T we have solution itself, and mapping from characteristic to real variables is more difficult

How do we solve hyperbolic problems?

Finite volume

Very closely related to finite differences. We look for a numerical method which mimics as much as possible the dynamics of the physical problem.

Elements in finite volume are called cells, which can be a □ or ▢ or anything else. In practice, we can have generic elements close to the boundary more structured mesh good grid around boundaries, caution away.

We'll stick to the 1D case:

Introduce: Usually called cell average. Want to find a good average value for the cell equation very similar to the real, original one. Generic numerical method related to finite volumes. Different methods depend on the choice of numerical flux. Most compute values on the boundaries can have discontinuities, values of U on two sides:

1. First problem: find a meaningful U
2. Second problem: missing time

Stability constraints! Choose quite often explicit method: but it is an unstable method. Must be then an upwind method derived from physical considerations.

Numerical flux for upwind method? Very good for implementation; no need to worry about its sign. Time loop & space loop (to update all values). To change the method, just change the numerical flux, not the structure of the code. To increase the accuracy of the method must increase the degree of polynomials in each cell.

Problem of discontinuity: shockwave controlling diffusion. Linear case initial solution transported. Non-linearity leads to break wave or reach a steady state. Non-linear case (e.g., Burger equation) can create discontinuity. Non-continuous initial conditions. Want to capture the shock (discontinuity) and propagate it correctly.

For simplicity, we'll consider 1D finite volumes in our analysis:

∂u/∂t + a ∂u/∂x = 0

u_iⁿ⁺¹ = u_iⁿ + λ (F_i+1/2^H - F_i-1/2^H), λ = Δt/Δx

Forward Euler Centered (FEC)

u_iⁿ⁺¹ - u_iⁿ + a (u_iⁿ - u_iⁿ) / Δx = 0

u_iⁿ⁺¹ - u_iⁿ - λ/2 a (u_iⁿ⁺¹ - u_iⁿ)

⇔ F_i+1/2^H = 1/2 [a (u_iⁱⁿ - u_i) ]

Note: Not the best one, because not respectful of the physics of the problem: would be more "example" to consider only upwind infos

Lax-Friedrichs (LF)

u_iⁿ⁺¹ = 1/2 (uⁿ⁺¹_i+1 + uⁿ_i-1) - λ/2 (uⁱⁿ_i - uⁱⁿ_i-1)

⇔ F_i+1/2^H = 1/2 [a (uⁱⁿ_i + uⁱⁿ_i-1) - λ (uⁱⁿ_i-1 - uⁱⁿ_i)]

Lax-Wendroff (LW)

• Compute u_iⁿ⁺¹ as Taylor-expanded in time: u_iⁿ⁺¹ = u_iⁿ + (∂u∕∂t)| _i Δt + (∂² u∕∂t²)| _i(Δt²/2 + (Δt³)

• From original equation: ∂u/∂t = -a ∂u/∂x ⇒ (∂²u∕∂t²) = -∂u∕∂t ∂x = -a ∂²u∕∂x²

uⁿ⁺¹ = uⁿ⁺¹ + (∂u/∂x) _i Δt + (a²∂²u/∂x²)| _i (Δt²/2 …

Discretizing with centered approximations:

u_iⁿ⁺¹ = u_iⁿ Δ(u_iⁿ - u_iⁱⁿa) Δt + a² Δ(u_iⁱⁿ -u_iⁱⁿa) (Δt²/2)

⇔ F_i+1/2^H = 1/2 [ a (uⁱⁿ_i + uⁱ_i) - λ² a²(uⁱⁿ_i-1 - uⁱⁿ_i)]

Upwind (up)

u_iⁿ⁺¹ - u_i^n-1 + a u_i^n-1 - u_iⁱⁿ - u_iⁿ)

• How do we choose which method to use? (In practice LW, UP are the most used; FEC almost fewer) (Check Consistency, stability, convergence)

• Consistency LTE (Local Truncation Error): "remainder of unexactness of the solution"

• FoL upwind

U(x,t) analytical solution T: T = (u(x_i, tⁿ⁺¹) - u(x_i, tⁿ)) + u(x_i, tⁿ) - u(x_i, t) )

• Function Error TF: τ = max _n,i | T|↳ If τ → ∞ when Δt , Δx → 0"

1 |If τ"(Δt, Δx) = O(Δt^p Δx^q) method is said to be of order p with respect to time and q with respect to space

Let us check consistency and order of upwind

In the same way, we can proceed to compute space approx, getting to O(Δx) q = 1

For upwind p = 1, q = 1

For FEC p = 1, q = 2

For LF O(Δx ₂ + Δt + o x₂)

For LW O(Δt ² + Δx ² + o x ₂Δt)

Convergence

lim max (u(x_i, tⁿ) - uⁿ) = 0

All of the methods we’ve seen are actually convergent, but most of the times it is better not to prove it directly but by means of: consistency + stability ⇒ convergence.

If method stable order of convergence is the same of the order of consistency

Stability

For T there is a C_T(T) such that ∀Δx > 0 there is δ(Δx) such that for 0c 0t c J₀, ‖u_n‖_Δ c C_T ‖u⁰‖_Δ for nQt ≤ T‖.‖_Δ is Norm delta; function inside is defined pointwise, so it is a vector containing [uⁿ, u_i…, uⁿ] all averages

Delta norm must be able to measure numerical solution. Problem with this is discretization of the norm: do not take into account all size reduction when number of cells increases

Δx is introduced to 'weight' value introduced in the computation

‖uⁿ‖_Δ discrete norm of discrete solution at time np = C_T + 1 method is said to be strongly stable ‖uⁿ‖_Δ ≤ ‖u^n-1‖_Δ true at each step; ‖uⁿ‖_Δ ≤ ‖u^n-1‖_Δ

Let us check stability of upwind method

G uⁿ⁺¹ = uⁿ + λΔ(uⁿ_i+1 - uⁿ_i) = (1- λΔ)

Numerical domain of dependence must contain Physical domain of dependence.

Courant-Friedrichs-Lewy (CFL) condition ➔ c = a Δt/Δx = Courant number

For systems, λ_i ∂Δt/Δx < 1

For FEC never strong stability, could have stability if Δt ∈ (Δx/σ)², with stability constant C = e^h/2

Another method to check stability is Von Neumann Method, using ‖uⁿ‖_l²∂uⁿ/∂t = a ∂u/∂x = 0

Numerical method will introduce some dissipation, so we’re ready to accept energy reduction

u_o(x_j) = Σ_k α_k e^ikx_j

Let's consider FEC

uⁿ⁺¹_j - uⁿ_j = a∂t/Δx (u^n_j+1 - u^n_j)

Direct amplification factor- y_jⁿ or yⁿ _j+1 can iterate, ending up with: yⁿ _j = (3v)ⁿ yⁿ _j

If I want the method to be stable, I need |U (n)|_l2≤ ||U(n-1)||_l2

Check how things work in this case:|r|^1/2 = (1 + Δt² Δx^-2 sin² (∇X))

Let's check for the upwind:

uⁿ⁺¹_j - uⁿ_j = Δt^νφνⁿ_xβ_w = αρξ^φc^SUPβ_w = 0 such that |r|^x,y

Numerical solution will obviously differ from numerical analysis in two ways:

1. Amplitude of the wave is reduced (*)
2. Peak is moved (forward or backward)

(*) - Phase modification

Upwind: function V that satisfies the equation ΔV + ΔV.x(t) - V_x(t)