GSoC 2016: Exponential Integrators Chiara Segala Mentor: Prof. Marco - PowerPoint PPT Presentation

GSoC 2016: Exponential Integrators Octave Conference 2017 March 21, 2017 GSoC 2016: Exponential Integrators Chiara Segala Mentor: Prof. Marco Caliari

GSoC 2016: Exponential Integrators Outline Exponential Integrators Matrix functions Implementation of solvers Numerical experiments

GSoC 2016: Exponential Integrators Exponential Integrators 1 Exponential Integrators 2 Matrix functions 3 Implementation of solvers 4 Numerical experiments

GSoC 2016: Exponential Integrators Exponential Integrators Exponential Integrators Parabolic problem u ′ ( t ) + Au ( t ) = f ( t ) , u ( t 0 ) = u 0

GSoC 2016: Exponential Integrators Exponential Integrators Exponential Integrators Parabolic problem u ′ ( t ) + Au ( t ) = f ( t ) , u ( t 0 ) = u 0 Solution (variation of constants formula) � h n u ( t n + 1 ) = e − h n A u ( t n ) + e − ( h n − τ ) A f ( t n + τ ) d τ 0

GSoC 2016: Exponential Integrators Exponential Integrators Exponential Integrators Parabolic problem u ′ ( t ) + Au ( t ) = f ( t ) , u ( t 0 ) = u 0 Solution (variation of constants formula) � h n u ( t n + 1 ) = e − h n A u ( t n ) + e − ( h n − τ ) A f ( t n + τ ) d τ 0 First numerical scheme (exponential quadrature rule) � s u n + 1 = e − h n A u n + h n b i ( − h n A ) f ( t n + c i h n ) i = 1 weights b i ( z ) can be expressed as linear combinations of the functions � 1 θ k − 1 e ( 1 − θ ) z ϕ k ( z ) = ( k − 1 )! d θ, k ≥ 1 . 0 ϕ k + 1 ( z ) = ϕ k ( z ) − ϕ k ( 0 ) ϕ k ( 0 ) = 1 ϕ 0 ( z ) = e z , k ! , z

GSoC 2016: Exponential Integrators Exponential Integrators Example s = 1 u n + 1 = e − h n A u n + h n ϕ 1 ( − h n A ) f ( t n + c 1 h n ) = u n + h n ϕ 1 ( − h n A )( f ( t n + c 1 h n ) − Au n )

GSoC 2016: Exponential Integrators Exponential Integrators Example s = 1 u n + 1 = e − h n A u n + h n ϕ 1 ( − h n A ) f ( t n + c 1 h n ) = u n + h n ϕ 1 ( − h n A )( f ( t n + c 1 h n ) − Au n ) Exponential Euler method for c 1 = 0 Butcher tableau 0 ϕ 1 Exponential midpoint rule for c 1 = 1 2

GSoC 2016: Exponential Integrators Exponential Integrators Exponential Runge-Kutta methods Semilinear problem u ′ ( t ) = F ( t , u ) = Au ( t ) + g ( t , u ( t )) , u ( t 0 ) = u 0

GSoC 2016: Exponential Integrators Exponential Integrators Exponential Runge-Kutta methods Semilinear problem u ′ ( t ) = F ( t , u ) = Au ( t ) + g ( t , u ( t )) , u ( t 0 ) = u 0 Numerical Exponential Runge-Kutta scheme � s u n + 1 = u n + h n b i ( h n A )( G ni + Au n ) i = 1 s � U ni = u n + h n a ij ( h n A )( G nj + Au n ) j = 1 G nj = g ( t n + c j h n , U nj )

GSoC 2016: Exponential Integrators Exponential Integrators Exponential Runge-Kutta methods Semilinear problem u ′ ( t ) = F ( t , u ) = Au ( t ) + g ( t , u ( t )) , u ( t 0 ) = u 0 Numerical Exponential Runge-Kutta scheme � s u n + 1 = u n + h n b i ( h n A )( G ni + Au n ) i = 1 s � U ni = u n + h n a ij ( h n A )( G nj + Au n ) j = 1 G nj = g ( t n + c j h n , U nj ) Reformulation s � u n + 1 = u n + h n ϕ 1 ( h n A ) F ( t n , u n ) + h n b i ( h n A ) D ni i = 2 i − 1 � U ni = u n + h n c i ϕ 1 ( c i h n A ) F ( t n , u n ) + h n a ij ( h n A ) D nj j = 2 D nj = g ( t n + c j h n , U nj ) − g ( t n , u n )

GSoC 2016: Exponential Integrators Exponential Integrators Exponential Rosenbrock methods u ′ ( t ) = F ( t , u ) , u ( t 0 ) = u 0

GSoC 2016: Exponential Integrators Exponential Integrators Exponential Rosenbrock methods u ′ ( t ) = F ( t , u ) , u ( t 0 ) = u 0 Numerical Exponential Rosenbrock scheme s � u n + 1 = u n + h n ϕ 1 ( h n J n ) F ( t n , u n ) + h 2 n ϕ 2 ( h n J n ) v n + h n b i ( h n J n ) D ni i = 2 i − 1 � U ni = u n + h n c i ϕ 1 ( c i h n J n ) F ( t n , u n ) + h 2 n c 2 i ϕ 2 ( c i h n J n ) v n + h n a ij ( h n J n ) D nj j = 2

GSoC 2016: Exponential Integrators Exponential Integrators Exponential Rosenbrock methods u ′ ( t ) = F ( t , u ) , u ( t 0 ) = u 0 Numerical Exponential Rosenbrock scheme s � u n + 1 = u n + h n ϕ 1 ( h n J n ) F ( t n , u n ) + h 2 n ϕ 2 ( h n J n ) v n + h n b i ( h n J n ) D ni i = 2 i − 1 � U ni = u n + h n c i ϕ 1 ( c i h n J n ) F ( t n , u n ) + h 2 n c 2 i ϕ 2 ( c i h n J n ) v n + h n a ij ( h n J n ) D nj j = 2 where J n = ∂ F ∂ u ( t n , u n ) v n = ∂ F ∂ t ( t n , u n ) g n ( t , u ) = F ( t , u ) − J n u − v n t D nj = g n ( t n + c j h n , U nj ) − g n ( t n , u n )

GSoC 2016: Exponential Integrators Matrix functions 1 Exponential Integrators 2 Matrix functions 3 Implementation of solvers 4 Numerical experiments

GSoC 2016: Exponential Integrators Matrix functions ϕ -functions � 1 θ k − 1 e ( 1 − θ ) z ϕ k ( z ) = ( k − 1 )! d θ, k ≥ 1 . 0 ϕ k + 1 ( z ) = ϕ k ( z ) − ϕ k ( 0 ) ϕ k ( 0 ) = 1 ϕ 0 ( z ) = e z , k ! , z

GSoC 2016: Exponential Integrators Matrix functions ϕ -functions � 1 θ k − 1 e ( 1 − θ ) z ϕ k ( z ) = ( k − 1 )! d θ, k ≥ 1 . 0 ϕ k + 1 ( z ) = ϕ k ( z ) − ϕ k ( 0 ) ϕ k ( 0 ) = 1 ϕ 0 ( z ) = e z , k ! , z ϕ -functions code: scaling and squaring + Padé approximation

GSoC 2016: Exponential Integrators Matrix functions Accuracy ϕ k ( ± 10 ) phikm( ± 10) phipade( ± 10,k) ϕ 1 ( 10 ) 2202.54657948068 2202.54657962449 ϕ 2 ( 10 ) 220.154657948068 220.154657962442 ϕ 3 ( 10 ) 21.9654657948068 21.9654657961828 ϕ 4 ( 10 ) 2.17987991281401 2.17987991294347 ϕ 1 ( − 10 ) 0.0999954600070237 0.0999954600070233 ϕ 2 ( − 10 ) 0.0900004539992977 0.0900004539992175 ϕ 3 ( − 10 ) 0.0409999546000702 0.0409999546000398 ϕ 4 ( − 10 ) 0.0125666712066596 0.0125666712066530

GSoC 2016: Exponential Integrators Matrix functions Augmented matrix Theorem (see Al-Mohy and Higham, 2011, Th. 2.1) Let A ∈ C n × n , W = [ w 1 , w 2 , . . . , w p ] ∈ C n × p , τ ∈ C , and � A � � 0 � W I p − 1 � ∈ C ( n + p ) × ( n + p ) , ∈ C p × p , A = J = 0 J 0 0 Then for X = ϕ l ( τ � A ) with l ≥ 0 we have j � τ k ϕ l + k ( τ A ) w j − k + 1 , X ( 1 : n , n + j ) = j = 1 : p . k = 1

GSoC 2016: Exponential Integrators Matrix functions Augmented matrix for exponential integrators W (: , p − k + 1 ) = u k , k = 1 : p , l = 0 , τ = t − t 0 � � e ( t − t 0 ) A A = X 12 A ) = e ( t − t 0 ) � X = ϕ 0 (( t − t 0 ) � e ( t − t 0 ) J 0 p � ϕ k (( t − t 0 ) A )( t − t 0 ) k u k X ( 1 : n , n + p ) = k = 1

GSoC 2016: Exponential Integrators Matrix functions Augmented matrix for exponential integrators W (: , p − k + 1 ) = u k , k = 1 : p , l = 0 , τ = t − t 0 � � e ( t − t 0 ) A A = X 12 A ) = e ( t − t 0 ) � X = ϕ 0 (( t − t 0 ) � e ( t − t 0 ) J 0 p � ϕ k (( t − t 0 ) A )( t − t 0 ) k u k X ( 1 : n , n + p ) = k = 1 p � u ( t ) = e ( t − t 0 ) A u 0 + ϕ k (( t − t 0 ) A )( t − t 0 ) k u k ˆ k = 1 = e ( t − t 0 ) A u 0 + X ( 1 : n , n + p ) A � u 0 � 0 ] e ( t − t 0 ) � = [ I n e p

GSoC 2016: Exponential Integrators Matrix functions Augmented matrix for exponential integrators W (: , p − k + 1 ) = u k , k = 1 : p , l = 0 , τ = t − t 0 � � e ( t − t 0 ) A A = X 12 A ) = e ( t − t 0 ) � X = ϕ 0 (( t − t 0 ) � e ( t − t 0 ) J 0 p � ϕ k (( t − t 0 ) A )( t − t 0 ) k u k X ( 1 : n , n + p ) = k = 1 p � u ( t ) = e ( t − t 0 ) A u 0 + ϕ k (( t − t 0 ) A )( t − t 0 ) k u k ˆ k = 1 = e ( t − t 0 ) A u 0 + X ( 1 : n , n + p ) A � u 0 � 0 ] e ( t − t 0 ) � = [ I n e p � � A � �� � u 0 η W η = 2 −⌈ log 2 ( � W � 1 ) ⌉ u ( t ) = [ I n ˆ 0 ] exp ( t − t 0 ) , η − 1 e p 0 J

GSoC 2016: Exponential Integrators Matrix functions Double augmented matrix Theorem (Double augmented matrix) Let A ∈ C n × n , W = [ w 1 , w 2 , . . . , w p ] ∈ C n × p , V = [ v 1 , v 2 , . . . , v q ] ∈ C n × q , τ ∈ C , and � A � W V � � ∈ C ( n + p + q ) × ( n + p + q ) , J ∈ C p × p , K ∈ C q × q . A = 0 J 0 0 0 K Then for X = ϕ l ( τ � � A ) with l ≥ 0 we have j � τ k ϕ l + k ( τ A ) w j − k + 1 , X ( 1 : n , n + j ) = j = 1 : p . k = 1 i − p � τ k ϕ l + k ( τ A ) v i − k + 1 − p , X ( 1 : n , n + i ) = i = ( p + 1 ) : ( p + q ) . k = 1

GSoC 2016: Exponential Integrators Matrix functions expmv e A B ≈ ( T m ( s − 1 A )) s B m = ⇒ select_taylor_degree

GSoC 2016: Exponential Integrators Matrix functions expmv e A B ≈ ( T m ( s − 1 A )) s B m = ⇒ select_taylor_degree function [f] = expmv (A, b) . . . M = select_taylor_degree (A, b); . . .

GSoC 2016: Exponential Integrators Matrix functions expmv e A B ≈ ( T m ( s − 1 A )) s B m = ⇒ select_taylor_degree function [f] = expmv (A, b) . . . M = select_taylor_degree (A, b); . . . M = select_taylor_degree (A, b); f = expmv (A, b, M);

GSoC 2016: Exponential Integrators Matrix functions expmv e A B ≈ ( T m ( s − 1 A )) s B m = ⇒ select_taylor_degree function [f] = expmv (A, b) . . . M = select_taylor_degree (A, b); . . . M = select_taylor_degree (A, b); f = expmv (A, b, M); ⇒ � A k � 1 select_taylor_degree =