Transcript Monte Carlo for Solving Linear Operator Equations

Monte Carlo for PDEs
Fall 2013
Review
• Last Class
– Monte Carlo Linear Solver
• von Neumann and Ulam method
• Randomize Stationary iterative methods
• Variations of Monte Carlo solver
• This Class
– Fredholm integral equations of the second kind
– The Dirichlet Problem
– Eigenvalue Problems
• Next Class
– Green’s Function Monte Carlo
Fredholm integral equations of the
second kind
• The integral equation
𝑓 𝑥 =𝑔 𝑥 +
𝐾 𝑥, 𝑦 𝑓 𝑦 𝑑𝑦
may be solved by Monte Carlo methods.
Since the integral can be approximated by a quadrature formula:
𝑁
𝑏
𝑦 𝑥 𝑑𝑥 =
𝑎
𝑤𝑗 𝑦 𝑥𝑗
𝑖=1
Fredholm integral equations of the
second kind
• The integral equation can be transformed to be
𝑁
𝑓 𝑥 =𝑔 𝑥 +
𝑤𝑗 𝐾 𝑥, 𝑦𝑗 𝑓 𝑦𝑗
𝑗=1
evaluate it at the quadrature points:
𝑁
𝑓 𝑥𝑖 = 𝑔 𝑥𝑖 +
𝑤𝑗 𝐾 𝑥𝑖 , 𝑦𝑗 𝑓 𝑦𝑗
𝑗=1
Let 𝑓 be the vector 𝑓 𝑥𝑖 , 𝑔 be the vector 𝑔 𝑥𝑖 and 𝐾 be the
matrix 𝑤𝑗 𝐾 𝑥𝑖 , 𝑦𝑗 , the integral equation becomes
𝑓 = 𝐾𝑓 + 𝑔
where 𝑓 is the unknown vector.
The Dirichlet Problem
• Dirichlet’s problem is to find a function 𝑢, which is continuous and
differentiable over a closed domain 𝐷 with boundary 𝐶, satisfying
𝛻 2 𝑢 = 0 𝑜𝑛 𝐷,
𝑢 = 𝑓 𝑜𝑛 𝐶.
where 𝑓 is a prescribed function, and
operator.
𝛻2
=
𝜕2 𝑢
𝜕𝑥
+
𝜕2 𝑢
𝜕𝑥
is the Laplacian
Replacing 𝛻 2 by its finite-difference approximation,
1
𝑢 𝑥, 𝑦 = 𝑢 𝑥, 𝑦 + ℎ + 𝑢 𝑥, 𝑦 − ℎ + 𝑢 𝑥 + ℎ, 𝑦 + 𝑢 𝑥 − ℎ, 𝑦
4
The Dirichlet Problem
• Suppose the boundary 𝐶 lies on the mesh, the previous
equations can be transformed into
𝑢 = 𝐻𝑢 + 𝑓
– The order of 𝐻 is equal to the number of mesh points in 𝐷.
1
– 𝐻 has four elements equal to in each row corresponding to an
4
interior point of 𝐷, all other elements being zero.
– 𝑓 has boundary values corresponding to an boundary point of 𝐶, all
other interior elements being zero.
– The random walk starting from an interior point 𝑃, terminates when it
hits a boundary point 𝑄. The 𝑓(𝑃) is an unbiased estimator of 𝑢(𝑄).
Eigenvalue Problems
• For a given 𝑛 × 𝑛 symmetric matrix 𝐻
𝐻𝑥𝑖 = 𝜆𝑖 𝑥𝑖 , 𝑥𝑖 ≠ 0
assume that 𝜆1 > 𝜆2 ≥ ⋯ ≥ 𝜆𝑛 , so that 𝜆1 is the dominant
eigenvalue and 𝑥1 is the corresponding eigenvector.
For any nonzero vector 𝑢 = 𝑎1 𝑥1 + 𝑎2 𝑥2 + ⋯ + 𝑎𝑛 𝑥𝑛 , according
to the power method,
𝐻𝑘 𝑢
lim
= 𝑎1 𝑥1
𝑘→∞ 𝜆 𝑘
1
We can obtain a good approximation of the dominant
eigenvector of 𝐻 from the above.
Eigenvalue Problems
Similar to the idea behinds Monte Carlo solver that
𝐻 𝑘 𝑢 = 𝑖1 … 𝑖𝑘 ℎ 𝑖0𝑖1 ℎ 𝑖1𝑖2 … ℎ 𝑖𝑘−1𝑖𝑘 a 𝑖𝑘
=
𝑖1 …
𝑖𝑘 𝑝 𝑖0 𝑖1 𝑝 𝑖1 𝑖2
…𝑝 𝑖
𝑘−1 𝑖𝑘
𝑝𝑖 𝑣
𝑘
𝑖0 𝑖1
𝑣 𝑖1𝑖2 … 𝑣 𝑖𝑘−1𝑖𝑘 u 𝑖𝑘 /p 𝑖
𝑘
we can do sampling on 𝐻 𝑘 𝑢 to estimate its value, and then
evaluate the dominant eigenvector 𝑥1 by a proper scaling.
From the Rayleigh quotient,
𝑥 𝑇 𝐻𝑥
𝜆= 𝑇
𝑥 𝑥
the dominant eigenvalue 𝜆1 be approximated based on the
estimated vector of 𝑥1 .
Summary
• This Class
– Monte Carlo Linear Solver
• von Neumann and Ulam method
• Randomize Stationary iterative methods
• Variations of Monte Carlo solver
– Fredholm integral equations of the second kind
– The Dirichlet Problem
– Eigenvalue Problems
What I want you to do?
• Review Slides
• Work on Assignment 4