A Monte Carlo method to compute principal eigenelements of some - PDF document

A Monte Carlo method to compute principal eigenelements of some linear operators Organization of the talk 1. Introduction 2. homogeneous neutron transport operators 3. Inhomogeneous operators 4. Laplace operator 5. Variance reduction + eigenfunction

The Monte Carlo method Numerical computation of A = E ( X ) X random variable, X i samples of X Monte Carlo method: law of large numbers N A ≃ 1 � X i N i =1 Central-limit theorem: con�dence interval N N 1 X i − a σ X ≤ E ( X ) ≤ 1 X i + a σ X � � √ √ N N N N i =1 i =1 σ X Speed of convergence: √ N 1

Cauchy problem in a domain D ∂u ∂t = Lu Neutron transport: sign of the principal eigen- value α 1 determins whether the system is sub- critical or supercritical. Di�usions operators : speed of convergence towards the steady state. Numerical computation of α 1 by deterministic methods: discretization of L in many dimen- sions + power method. 2

Spectral expansion in L 1 ( D ) u ( x, v, t ) = D 0 e α 1 t + D 1 e α 2 t + ... + D k e α k t + o ( e α k t ) α 0 real and simple, D 0 independent of t Neutron transport: Mika (1968), Vidav (1971) Leading terms E ⊂ D 1 E u ( x, v, t ) dxdv ) ≃ 1 � t log( C 0 e α 1 t + C 1 e α 2 t ) t log( Stochastic representation 1 � lim t log( E u ( x, v, t ) dxdv ) = α 1 t →∞ Approximation model e ( α 2 − α 1 ) t 1 E u ( x, v, t ) dxdv ) ≃ α 1 +log( C 0 ) + C 1 � t log( . t C 0 t 3

Monte Carlo estimators of α 1 Approximation model e ( α 2 − α 1 ) t 1 E u ( x, t ) dx ) ≃ α 1 +log( C 0 ) + C 1 � t log( t C 0 t Interpolation method α 1 ( t 1 , t 2 ) = log( � E u ( x, t 2 ) dx ) − log( � E u ( x, t 1 ) dx ) t 2 − t 1 Linear least square problem q − 1 ( α 1 + β � E u ( x, v, t i ) dx )) 2 � log( t i t i i = p Optimization and variance reduction tools Maire,Talay IMAJNA(06), Lejay,Maire MCS(07) 4

Neutron transport equations Cauchy problem in A × V ⊂ ℜ 3 × ℜ 3 ∂u ∂t ( x, v, t ) = Tu ( x, v, t ) Neutron transport operator T Tu ( x, v, t ) = − v ∇ x u ( x, v, t ) − Σ( x, v ) u ( x, v, t ) + � V f ( x, z, v ) u ( x, z, t ) dz + s ( x, v, t ) u ( x, v, t ) density of neutrons Initial condition u ( x, v, 0) = u 0 ( x, v ) Boundary conditions Heterogeneous problem in dimension 7 5

Criticality problem Sources come from �ssion T λ u ( x, v, t ) = − v ∇ x u ( x, v, t ) − Σ( x, v ) u ( x, v, t ) + V f r ( x, z, v ) u ( x, z, t ) dz + 1 � � V f f ( x, z, v ) u ( x, z, t ) dz λ v. ∇ u : loss of neutrons by transport Σ u : loss of neutrons by collisions V f r udz : sources from slowing � 1 V f f udz : source from �ssion � λ Criticality : �nd λ such that the principal ei- genvalue α λ of the operator T λ is zero 6

Neutron transport operators Cauchy problem in A × V ∂u ( x, y, t ) = b ( x, y ) ∇ x u ( x, y, t ) + c ( x, y ) u ( x, y, t ) + ∂t �� � V u ( x, z, t ) Π x,y ( dz ) − u ( x, y, t ) γ ( x, y ) Feynman-Kac representation of u ( x, y, t ) �� t � �� E x,y 1 τ x,y A >t u 0 ( X t , Y t ) exp 0 c ( X s , Y s ) ds τ x,y exit time from A A Y t jump process on V X t is solution of dX t dt = b ( X t , Y t ) Exact simulation schemes 7

The model of Lehner and Wing � 1 � � ∂u ∂t = − v∂u 1 − 1 u ( x, v ′ , t ) dv ′ − u ( x, v, t ) ∂x + c 2 + ( c − 1) u ( x, v, t ) Stochastic representation of u ( x, v, t ) �� t � � � u 0 ( X t , Y t ) exp 0 ( c − 1) ds 1 t<σ x,v E x,v D σ x,v : exit time from D = [0 , A ] D dX t dt = − Y t ; X 0 = x and Y 0 = v . Velocity at collision uniform law on [ − 1 , 1] . Time between two collisions � t � � F x,v ( t ) = 1 − exp = 1 − exp ( − ct ) − 0 cds 8

Choice of u 0 : u 0 ( x, v ) ≡ 1 u ( x, v, t ) = P x,v ( σ x,v > t ) exp( c − 1) t D � � E P x,v ( σ x,v E u ( x, v, t ) dxdv = exp( c − 1) t > t ) dxdv D Monte-Carlo Approximation > t ) dxdv ≃ vol ( E ) N ( t ) � E P x,v ( σ x,v D N N ( t ) number of trajectories still alive at time t among N starting from a random uniform point in E. 1 E u ( x, v, t ) dxdv ) ≈ c − 1 + 1 t log( N ( t ) � t log( ) N 9

Naive approximation of α 0 0 ’r1c10mili’ -0.005 -0.01 -0.015 -0.02 -0.025 -0.03 -0.035 -0.04 0 50 100 150 200 250 300 350 400 t log( N ( t ) 1 N ) as a function of time t . 10 millions of simulations, A = 8 Variance increases as t increases Rare events simulations 10

Approximation of α 0 using interpolation c − 1 + 1 log( N ( t 1 ) ) ≃ α 0 + log( K 0 ) t 1 N t 1 c − 1 + 1 log( N ( t 2 ) ) ≃ α 0 + log( K 0 ) t 2 N t 2 α 0 ( t 1 , t 2 ) = c − 1 + log( N ( t 2 ) ) − log( N ( t 1 ) ) N N . t 2 − t 1 Con�dence interval for α 0 ( t 1 , t 2 ) 11

Approximation using interpolation -0.024 ’vp1c10mili20’ -0.026 -0.028 -0.03 -0.032 -0.034 -0.036 -0.038 -0.04 0 20 40 60 80 100 120 140 160 180 α 0 ( t 1 , t 1 + 20) − c + 1 as a function of t 1 . 10 millions of simulations Approximation stable and accurate 12

Criticality computation c = 1 . 036 binf bsup t 1 α 0 min ( t 1 , t 2 ) 50 -4.28 10 − 4 -4.12*10 − 4 -3.94 10 − 4 55 -4.28 10 − 4 -4.13*10 − 4 -3.94 10 − 4 60 -4.28 10 − 4 -4.16*10 − 4 -3.96 10 − 4 c = 1 . 037 binf bsup α 0 max ( t 1 , t 2 ) t 1 50 6.05 10 − 4 6.22 10 − 4 6.38 10 − 4 55 6.02 10 − 4 6.21 10 − 4 6.39 10 − 4 60 5.9 10 − 4 6.19 10 − 4 6.4 10 − 4 13

Secant method gives c c c c = cmin − αmin cmax − cmin αmax − αmin = 1 . 036399 Error : 10 − 5 on c c . Dahl-Sjostrand (1978). Least square approximation Find α 0 and β = log( K 0 ) minimizing q ( α 0 + β − 1 � E u ( x, v, t i ) dxdv )) 2 . � log( t i t i i = p Linear least square problem Improvement of the accuracy 14

Other inhomogeneous problems Anisotropic problem ∂u = − v∂u ∂t ∂x +( c − 1) u + c � V (1+3 µvv ′ ) u ( x, v ′ , t ) dv ′ . 2 ⇒ Change the law of the velocity 1 + 3 µvv ′ � � f ( v ′ ) = 2 Dimension 3 : Case of a sphere ∂u ( x, v, t ) = − v ∇ x u ( x, v, t ) + ( c − 1) u ( x, v, t ) + ∂t c ( 1 � u ( x, z 1 , t ) dz 1 − u ( x, z, t )) 4 π S 2 Same accuracy 15

Inhomogeneous cases Zone 1 and 3 : splitting zones Zone 2 and 4 : very absorbing zones Zone 5 : absorbing zone 16

Additional di�culties �� t � �� u ( x, y, t ) = E x,y 1 σ x,y A >t u 0 ( X t , Y t ) exp 0 c ( X s , Y s ) ds 1. Following the motion precisely 2. Computation of u at di�erents times 3. Handle change of zones 4. Initial condition u 0 chosen using another numerical method 5. Computation of u on the most splitting zone 17

Laplace operator with Dirichlet conditions Cauchy problem in D ∂u ( x, t ) = 1 2 △ u ∂t Feynman-Kac representation of u ( x, t ) � � E x 1 τ D >t u 0 ( B t ) If u 0 ( x ) ≡ 1 then u ( x, t ) = P x ( τ D > t ) D : exit time from D of B t τ x Simulation schemes: Euler Scheme with or without boundary correction, walk on spheres or rectangles 18

Homogeneous transport or Laplace operator Branching method to compute P x ( τ D > T ) Cutting time to reach into slices Markov property P x ( τ D > T ) = P x ( τ D > T 1 ) P π T 1 ( τ D > T − T 1 ) Particular approximation of π T 1 Law of X t assuming that τ D > t → < u 0 , ϕ ∗ E x [ u 0 ( X t ) /τ D > t ] = E x [ u 0 ( X t ); τ D > t ] 1 > − < 1 , ϕ ∗ E x [1; τ D > t ] 1 > ϕ ∗ for t large , we obtain 1 < 1 ,ϕ ∗ 1 > Variance reduction + eigenfunction 19

Conclusion Good accuracy on the principal eigenvalue on many test-cases. Method well adapted to high dimensions Sensitive to the simulation schemes Approximation of the principal eigenfunction in the case of the Laplace operator Stochastic representation of the principal eigen- value for some homogeneous transport opera- tor. Computation of the second leading eigenvalue? 20