Commit cbb635d8 by Christoph Ruegge

 %% Cell type:markdown id: tags: # Poisson equation in 2d ## Poisson equation Find $u \in C^2(\Omega)\cap C(\overline{\Omega})$ with $$\begin{cases} \Delta u = f & \text{in } \Omega \\ -\Delta u = f & \text{in } \Omega \\ u = 0 & \mbox{on } \partial \Omega \end{cases}$$ for $\Omega = (0,1)^2$ and $f \in C(\Omega)$. For $f = 0$, the equation is known as *Laplace equation*. ## Laplace operator $$\Delta u := \sum_{i=1}^d \frac{\partial ^2 u}{\partial x_i^2}$$ ### Discretization - Equidistant grid $\Omega_N = I_N \times I_N$, where $$I_N := \left\{ x^N_j \colon j = 0, \ldots, N-1 \right\} \subset (0, 1), \qquad x^N_j := \left( j + \frac{1}{2} \right) h_N, \qquad h_N := \frac{1}{N}.$$ - Approximation of $\frac{\partial^2 u}{\partial x^2}$: $$\frac{u(x-h, y) - 2 u(x, y) + u(x+h, y)}{h^2} = \frac{\partial^2 u}{\partial x^2} (x, y) + \mathcal{O}(h^2),$$ and similarly for $\frac{\partial^2 u}{\partial y^2}$. - Approximation of Laplace operator: $$\frac{u(x, y-h) + u(x-h, y) - 4 u(x, y) + u(x, y+h) + u(x+h, y)}{h^2} = \Delta u(x, y) + \mathcal{O}(h^2)$$ - Define $$u^N_{jk} := u(x^N_j, x^N_k), \qquad f^N_{jk} := f(x^N_j, x^N_k),$$ for $j, k = 0, 1, \ldots, N-1$, and $$(\Delta^N u^N)_{jk} := \frac{u^N_{j,k-1} + u^N_{j-1,k} - 4 u^N_{jk} + u^N_{j,k+1} + u^N_{j+1,k}}{h_N^2}.$$ Here, $u^N_{-1,k} = u^N_{N,k} = u^N_{j,-1} = u^N_{j,N} = 0$ (boundary condition). - Solve $$\Delta^N u^N = f^N.$$ $$-\Delta^N u^N = f^N.$$ ## Code For performance reasons, the code is not written in the notebook. See files 01-hmatrix.hs, 02-repa.hs and 03-stencils.hs. ... ...
 ... ... @@ -8,7 +8,7 @@ import Numeric.LinearAlgebra.Helpers import Numeric.LinearAlgebra.Repa (linearSolve) poisson :: Int -> (Double -> Double -> Double) -> Mat Double poisson n f = computeS . reshape gridshape . fromJust . linearSolve laplace $asColumn rhs poisson n f = computeS . reshape gridshape . fromJust . linearSolve negLaplace$ asColumn rhs where -- coordinates h = 1 / fromIntegral n :: Double ... ... @@ -34,8 +34,8 @@ poisson n f = computeS . reshape gridshape . fromJust . linearSolve laplace $as -- reshape gridshape :: Array r DIM1 Double -> Array D DIM2 Double matshape = Z:.(n*n):.(n*n) :: DIM2 laplace :: Mat Double laplace = computeS$ fromFunction matshape computeElem negLaplace :: Mat Double negLaplace = computeS $fromFunction matshape computeElem where computeElem :: DIM2 -> Double computeElem (Z:.(!j):.(!k)) | j == k = 4*c ... ...  ... ... @@ -5,7 +5,7 @@ import Data.Array.Repa import Data.Array.Repa.Helpers (runCG, saveArrayAsPNG) poisson :: Int -> (Double -> Double -> Double) -> IO (Array U DIM2 Double) poisson n f = runCG laplace =<< computeP rhs poisson n f = runCG 1e-6 negLaplace =<< computeP rhs where -- coordinates h = 1 / fromIntegral n :: Double ... ... @@ -16,8 +16,8 @@ poisson n f = runCG laplace =<< computeP rhs rhs = fromFunction gridshape$ \(Z:.(!jx):.(!jy)) -> f (unsafeLinearIndex xs jx) (unsafeLinearIndex xs jy) -- laplace operator laplace :: Array U DIM2 Double -> Array D DIM2 Double laplace arr = fromFunction ext computeElem negLaplace :: Array U DIM2 Double -> Array D DIM2 Double negLaplace arr = fromFunction ext computeElem where !ext = extent arr get !idx | inShape ext idx = unsafeIndex arr idx ... ...
 ... ... @@ -35,7 +35,7 @@ import Data.Array.Repa.Stencil import Data.Array.Repa.Stencil.Dim2 poisson :: Int -> (Double -> Double -> Double) -> IO (Array U DIM2 Double) poisson n f = runCG laplace =<< computeP rhs poisson n f = runCG 1e-6 negLaplace =<< computeP rhs where -- coordinates h = 1 / fromIntegral n :: Double ... ... @@ -46,7 +46,7 @@ poisson n f = runCG laplace =<< computeP rhs rhs = fromFunction gridshape $\(Z:.(!jx):.(!jy)) -> f (unsafeLinearIndex xs jx) (unsafeLinearIndex xs jy) -- laplace stencil laplaceStencil = makeStencil2 3 3 computeElem :: Stencil DIM2 Double negLaplaceStencil = makeStencil2 3 3 computeElem :: Stencil DIM2 Double where computeElem (Z:. -1 :. 0) = Just$ -c computeElem (Z:. 1 :. 0) = Just $-c computeElem (Z:. 0 :. -1) = Just$ -c ... ... @@ -55,8 +55,8 @@ poisson n f = runCG laplace =<< computeP rhs computeElem _ = Nothing c = 1/h**2 -- laplace operator laplace :: Array U DIM2 Double -> Array D DIM2 Double laplace = delay . mapStencil2 (BoundConst 0) laplaceStencil negLaplace :: Array U DIM2 Double -> Array D DIM2 Double negLaplace = delay . mapStencil2 (BoundConst 0) negLaplaceStencil -- There is special syntax for stencils. The laplaceStencil above can be -- produced more easily by ... ...