Hyperbolic Problems

NUMERICAL METHODS FOR DIFFERENTIAL EQNS

HYPERBOLIC PDEs

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.

⇒ du/dt + a du/dx = 0

f(u) = 1/2 u² ⇒ Burgers

_∂u/_∂t + 2u _∂u/_∂x

a(u) = x0 ⇒ dx/dt

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

⇒ Sub = x₀

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

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

➔ Inflow is where a·n < 0

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

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

NUMERICAL METHODS FOR DIFFERENTIAL EQNS

HYPERBOLIC PDEs

Totally different nature with respect to others.

Let's start our analysis from conservation laws:

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

b c

a c₂

b c

general could be f(u,x)

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

NON CONSERVATIVE FORM

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

usually = 0

sink (0)

• if a(u) = const. ⇒ ∂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

∂x/∂t = a(u) ⇒ 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), i∈ℝ

Condition for system to be hyperbolic is A can be reduced to diag through T.

^T_i^-1 ∂^{v/ _∂t} + ∇ΛΓ^-1 ∂^{v/ _∂x} < 0 -> 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

1λ >0

2λ >0

3λ <0

->

b.c. 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 dynamics of the physic problem
• elements in finite volume are called cells □ 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 1D case:

Introduce:

Usually called cell average

Want to find good average value for the cell

eqn 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 2 sides:

1st problem: find a meaningful U

2nd problem: missing time

STABILITY CONSTRAINTS!

choose quite often explicit method:

but it is an unstable method

much 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 method just change numerical flux, not structure of the code

To increase accuracy of the method must increase 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 eqn) can create discontinuity.

non-continuous i.c.

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 eqn: ∂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 wrt time and q wrt space

Let us check consistency and order of upwind

In 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: consist. + stabil. ⇒ converg.

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 n

p = 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- λΔ) c

Num. dom. of dep. must contain Phys. dom. of dep.

Courant Friedrich Levy (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)

dir. 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 solutionwill obviously differ fromnumerical analysis in two ways:

1. 1) Amplitude ofthe wave is reduced (*)
2. 2) Peak is moved(forward or backward)(*) - Phase modification

- Upwind:

function V that satisfies equation

ΔV + ΔV.x(t) - V_x(t)