The Fast Sweeping Method for Convex Hamilton-Jacobi Equations and - PowerPoint PPT Presentation

The Fast Sweeping Method for Convex Hamilton-Jacobi Equations and Beyond Hongkai Zhao UC Irvine

Highlights of the fast sweeping method (FSM) ◮ Simplicity: a Gauss-Seidel type of iterative method. ◮ Optimal complexity: finite number of iterations independent of mesh size. ◮ A truly nonlinear method. ◮ The idea and iterative approach can be extended to more general problems.

A simple Hamilton-Jacobi (HJ) equation: the eikonal equation Isotropic case: �∇ u ( x ) � = c ( x ) , u ( x ) = 0 , x ∈ S The viscosity solution u ( x ) is the first arrival time for a front starting at 1 S with a propagation speed v ( x ) = c ( x ) . T=0 v n =f T 1 T 2 1 2 = 1 , Anisotropic case: [ ∇ u ( x ) M ( x ) ∇ u ( x )] M ( x ) is a SPD matrix HJ equation is nonlinear hyperbolic PDE. The proper definition of weak solution is called the viscosity solution.

Different interpretation of viscosity solution (1) u is a viscosity subsolution ( supersolution ) if for all φ ∈ C ∞ 0 (Ω) with u − φ attaining a local maximum (minimum) at some x 0 ∈ Ω then H ( x 0 , ∇ φ ( x 0 )) ≤ 0 ( ≥ 0) A viscosity solution is both a viscosity sub- and supersolution. φ u u φ (a) subsolution (b) supersolution ◮ A general definition using test function. Useful for mathematical proof. ◮ Non-constructive and difficult to use for designing numerical schemes.

Different interpretation of viscosity solution (2) Let u ǫ be the solution to H ( x , ∇ u ǫ ( x )) = ǫ ∆ u ǫ ( x ) , u ǫ → u as ǫ → 0. u is the viscosity solution. ◮ Singular perturbation ⇒ equation type is changed ⇒ more global dependence. ◮ Using numerical viscosity ⇒ Lax-Friedrichs scheme. ◮ General and useful for non-convex Hamiltonian.

Different interpretation of viscosity solution (3) Convex Hamiltonian: the viscosity solution is the optimal control solution u ( x ) = inf y ∈ ∂ Ω ( g ( y ) + l ( x , y )) � 1 � � 0 ρ ( ξ ( t ) , − ξ ′ ( t )) dt | ξ ∈ C 0 , 1 ([0 , 1] , Ω) , ξ (0)= x , ξ (1)= y l ( x , y )=inf ξ ( t ) where ρ ( x , q ) = max { p | H ( x , p )=0 } < p , q > . Ω Γ f x α ⋅ ( , ( )) x ◮ Characteristics for the hyperbolic PDE ⇒ local optimal paths. ◮ When characteristics intersect, take the optimal value among all characteristics ⇔ viscosity solution. ◮ If Ω is convex and H is homogeneous in space ⇒ the optimal path is a straight line and l ( x , y ) = max { p | H ( p )=0 } ( x − y ) · p .

Numerical methods Two crucial ingredients: ◮ Appropriate discretization (local solver). ◮ converges to the right viscosity solution (monotone), ◮ follows characteristics (upwind), ◮ deals with non-smoothness and has high order accuracy if possible (nonlinear). ◮ Solve a large nonlinear system by iterative methods: ◮ Jacobi iteration (pseudo time marching): explicit but need many iterations due to the finite speed of propagation and the CFL condition. ◮ Gauss-Seidel iteration: efficient nonlinear solver based on upwind monotone schemes and causality of the PDE.

Other related works ◮ Danielson’s algorithm (1980). ◮ Schneider, et al, post sweeping (1992). ◮ Bou´ e and Dupuis, fast sweeping method for stochastic control (1999). ◮ Jameson and Caughey use alternating sweeping to solve linearized steady compressible Euler equation (2001). ◮ Tsitsiklis (1995), Sethian (1996) fast marching method. Sethian and Vladimirsky (2001) ordered upwind scheme.

The FSM on rectangular grids Solve |∇ u ( x ) | = c ( x ), with u ( x ) = given , x ∈ S . ◮ Local solver ( nonlinear equation ): upwind and causal [( u i , j − u xmin ) + ] 2 +[( u i , j − u ymin ) + ] 2 = c 2 i , j , i , j = 1 , 2 , . . . (1) h h u xmin =min( u i − 1 , j , u i + 1 , j ) , u ymin =min( u i , j − 1 , u i , j + 1 ), causality, ◮ Initial guess: u i , j = u ( x i , j ) = given , x i , j ∈ S , u i , j = ∞ , x i , j / ∈ S . ◮ G-S iterations with four alternating orderings: (1) i = 1 : I , j = 1 : J (2) i = I : 1 , j = 1 : J (3) i = I : 1 , j = J : 1 (4) i = 1 : I , j = J : 1 Let u be the solution to (1) using current values of u i ± 1 , j , u i , j ± 1 , u new = u = min( u old i , j , u ) , control interpretation. i , j

The FSM in 1D In 1D, there are two directions of characteristics, left or right. Two sweeps are needed to compute the numerical solution. the distance function to two points x 1 , x 2 : | u x | = 1 , u ( x 1 ) = u ( x 2 ) = 0. the computed solution after first sweeping: left to right the computed solution after second sweeping: right to left

The fast sweeping algorithm in 2D Facts: ◮ The characteristics have infinitely many directions. ◮ However all directions can be classified into four groups: up-right, up-left, down-left and down-right. Each sweeping order covers one group of characteristics. i=1:I, j=1:J i=I:1, j=1:J i=1:I, j=1:J i=I:1, j=1:J i=1:i, j=J:1 i=I:1, j=J:1 j j i=1:I, j=1,J i=I:1, j=1:J i=I:1, j=J:1 i=1:I, j=J:1 i=I:1, j=J:1 i=1:i, j=J:1 i i (a) distance to a point (b) distance to a circle

Convergence and error estimate Theorem 1: The solution of the fast sweeping algorithm converges monotonically to the solution of the discretized system. d ( x , S ): the distance function to S .; u h ( x , S ): the numerical solution. Theorem 2: For a single data point S = { x 0 } , u h ( x , x 0 ) converges after 2 n sweeps in n dimensions and d ( x , x 0 ) ≤ u h ( x , x 0 ) ≤ d ( x , x 0 ) + O ( h | log h | ) (sharp!) Theorem 3: For arbitrary S in n dimensions, after 2 n iterations. u ( x i , j , S ) ≤ u h ( x i , j , S ) ≤ d ( x i , j , S ) + O ( | h log h | ) , where u ( x i , j , S ) is the solution to the discretized system.

General eikonal equation |∇ u ( x ) | = c ( x ), denote H ( p , x ) = | p | − c ( x ), where p = ∇ u . The characteristic equations are: x = ∇ p H = ∇ u  ˙ c ( x )  p = −∇ x H = ∇ c ( x ) ˙ u = ∇ u · ˙  ˙ x = c ( x ) Each characteristic can be segmented into a finite number of pieces and can be covered by the G-S sweeps successively. ! 5 1 2 x 0 " 3 4

Total direction variation along a characteristic  � ∇ u � � � | ˙ x | = � = 1 ,   � �  c ( x )  �  x = ∇ ˙ c ( x ) − ∇ u u · ∇ c x = ( I − P n ) ∇ c ( x )   ¨ · ˙ c ( x ) .    c c ∇ u P n is the projection on n = |∇ u | . ¨ x is the curvature. |∇ c ( x ) | Let K = max x ∈ Ω c ( x ) , D = diameter of the domain, and c M ( c m ) is the maximum (minimum) of c ( x ). The total direction variation is bounded by � � |∇ c ( x ) | � ds ≤ K length(Γ) ≤ DK c M | ¨ x | ds ≤ ds ≤ K c ( x ) c m Γ Γ Γ ⇒ The maximum number of turns is bounded by DK c M c m . 2 π

Maximum length of a characteristic Let the characteristics Γ joins x 0 ∈ S and a x ∈ Ω. � x � � c m ds ≤ c ( s ) ds = u ( x ) ≤ c ( s ) ds ≤ c M � x − x 0 � Γ Γ x 0 Hence � ds ≤ Dc M length(Γ) = c m Γ where D is the radius of the domain and c M ( c m ) is the maximum (minimum) of c ( x ). The maximum number of turns is bounded by DK c M c m . 2 π x x 0 � S

FSM on unstructured mesh ( Luo, Qian, Zhang & Zhao) Two key ingredients: ◮ Montone upwind scheme. ◮ Systematic ordering covering all characteristics.

Ordering on unstructured mesh The Key point: Systematic ordering such that each characteristics can be covered in a finite number of sweep orderings. ◮ Order all vertices according to the distance to a few reference points and sweep back and forth. characteristic reference point C reference point reference point A reference point B characteristic

Discretization using optimal control (Bornemann & Rasch) Piecewise linear approximation on a local mesh, Ω h C , around the vertex C . ◮ Local constant approximation of the Hamiltonian: H h C ( p ) = H ( C , p ). ◮ Local optimal control problem: u h ( C ) = � l h C ( C , y ) + g h � inf C ( y ) . y ∈ ∂ Ω h C C ( y ) is the piecewise linear interpolation of u h at vertices on ∂ Ω h ◮ g h C . ◮ The optimal path is a straight line! l h C ( x , y ) = max C ( p )=0 } ( x − y ) · p . { p | H h C A B Mesh Ω h around vertex C

Properties of the discretization The discretization is of the implicit form: F h C ( C , u h C , u h A , u h B ) = 0. The scheme is ◮ consistent H h if u is linear, ∇ u = p , C ( C , ∇ u ) = H ( C , p ) , ◮ upwind and monotone 0 ≤ ∂ u h , ∂ u h ∂ u h + ∂ u h C C C C ≤ 1 , = 1 ∂ u h ∂ u h ∂ u h ∂ u h A B A B ◮ For anisotropic eikonal equation: [ ∇ u ( x ) M ( x ) ∇ u ( x )] 1 / 2 = 1, the implicit relation can be explicitly solved (Kao, Osher & Tsai 05, Bornemann & Rasch 06, Qian, Zhang & Z. 07) u h C = u h C ( C , u h A , u h B ) .

Convergence Theorem ◮ There exists a unique solution to the discretized nonlinear system. ◮ The fast sweeping iteration converges monotonically to the unique solution. ◮ The unique solution to the discretized solution converges to the viscosity solution as h → 0.

Two special cases on rectangular grids. C C (a) Five point stencils. (b) Nine point stencils. Nine point stencils will yield more accurate solutions due to better angle resolution.