Commit f66f1a08 authored by Jochen Schulz's avatar Jochen Schulz
Browse files

updated todo

parent a150cf87
# vim:ft=todo
LECTURE:
-exercise 23 : stiff nonstiff
-exercise 24: fft tiefpassfilter, jacobi mit repa, mandelbrot-irgendwas, poisson 1d
23: implicit euler in runge kutta einfuehren
stiff nonstiff
-exercise 23 fehlt
-exercise 24 fehlt
vorlesung 25: stencil muss auch erklaert werden (anhand differenzenquotient)
1d poisson beispiel fuer cg-loeser
25 poisson equation (2d statisch) mit repa (stencil!)
exercise 25: poisson 2d repa
25: stencil muss auch erklaert werden
25 poisson equation (2d statisch) mit/ohne repa (stencil!)
27 continous optimization (minimierung: russel? siam -probleme )
26 continous optimization (minimierung: russel? siam -probleme )
projektionsverfahren (russell) sudoku-solver (bfgs fertig machen)
exercise 26
exercise inverse
xx diskrete optimierung (lintim) shortest path implementation !
uebungen fuer die beispiele erstellen (teilweise teile in der vorlesung auslassen damit sie in der uebung gemacht werden koennen)
......
%% Cell type:code id: tags:
%% Cell type:markdown id: tags:
``` haskell
# poisson equation in 2d
```
%% Cell type:markdown id: tags:
- Poisson Problem also describes the stationary heat equation.
**Poisson equation** <br />
Find $u \in C^2(\Omega)\cap C(\overline{\Omega})$ with
$$\left \{ \begin{array}{rcll}- \triangle u & = & f & \mbox{in } \Omega\\
u & = & 0 & \mbox{auf } \partial \Omega\\ \end{array} \right.$$
for $\Omega=(0,1)^2$ and $f \in C(\Omega)$.
With $f = 0$ the equation is known as *laplace equation*
**Laplacian** <br />
$\triangle u := \sum_{i=1}^d \frac{\partial ^2 u}{\partial x_i^2}$
- equidistant gridsize: $h= \frac 1 N$,
$N \in \mathbb{N}$
- set of all gridpoints
$$Z_h := \left\{ (x,y) \in \overline{\Omega} \ \mid \ x=z_1h, \ y=z_2h \text{ mit } z_1,z_2 \in \mathbb{Z} \right\}.$$
- inner gridpoints: $\omega_h := Z_h \cap \Omega$
- approximation of $\frac{\partial ^2 u}{\partial x^2} (x,y)$
$$\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)$$
- approximation of $\frac{\partial ^2 u}{\partial y^2} (x,y)$
$$\frac{u(x ,y-h) - 2 u(x,y) + u(x,y+h)}{h^2} = \frac{\partial ^2 u}{\partial
y^2} (x,y) + \mathcal{O}(h^2)$$
- the sum gives the approximation for $\triangle u(x,y)$:
$$\frac{1}{h^2} \left( u(x,y-h) + u(x-h,y) - 4 u(x,y) + u(x,y+h) +
u(x+h,y) \right)$$
- definition $u_{i,j}:=u(ih,jh)$ ergibt an Gitterpunkten $(ih,jh)$
$$-u_{i,j-1} - u_{i-1,j} + 4 u_{i,j} - u_{i+1,j} - u_{i,j+1} = h^2 f_{ij}$$
with $i,j \in \{ 1, \dots , N-1 \}$ and $f_{ij}:=f(ih,jh)$.
- boundary conditions
$u_{0,i}=u_{N,i}=u_{i,0}=u_{i,N}=0$, $i=0, \dots ,N$.
- save 2D-values of $u$ in a vector: shape the inner variables
$$\begin{array}{cccc}
u(h,(N-1)h) & u(2h,(N-1)h) & \ldots & u((N-1)h,(N-1)h)\\
\vdots & \vdots & \vdots & \vdots \\
u(h,2h) & u(2h,2h) & \ldots & u((N-1)h,2h)\\
u(h,h), & u(2h,h) & \ldots & u((N-1)h,h)\\
\end{array}$$
this gives the vector $U_{i+(N-1)(j-1)}=u_{i,j}$.
Linear system for $U=(U_i)_{i=1}^{(N-1)^2}$
$$A U = F$$
with
- $F:=(f_i)_{i=1}^{(N-1)^2}$ mit $f_{i+(N-1)(j-1)}=f(ih,jh)$, $i,j \in \{1,
\dots ,N-1 \}$,
- $$A := \frac{1}{h^2} tridiag(-I_{N-1}, T, -I_{N-1}) \in \mathbb{R}^{(N-1)^2
\times (N-1)^2},$$
$$T := tridiag(-1,4,-1) \in \mathbb{R}^{(N-1)\times (N-1)}.$$
%% Cell type:code id: tags:
``` haskell
:extension NoMonomorphismRestriction
import Numeric.Container
import Numeric.LinearAlgebra.Algorithms
import Numeric.LinearAlgebra.Util
import Control.DeepSeq (force)
idx2grid :: Int -> Int -> (Int, Int)
idx2grid n i = (i `mod` n, i `div` n)
idxf :: Int -> (Int, Int) -> Double
idxf n (i, j)
| i == j = 4
| xi == xj && abs(yi - yj) == 1 = -1
| yi == yj && abs(xi - xj) == 1 = -1
| otherwise = 0
where (xi, yi) = idx2grid n i
(xj, yj) = idx2grid n j
poisson :: ([(Double, Double)], Matrix Double, Vector Double)
poisson = (mesh, reshape (n-1) $ head . toColumns $ loes, f)
where
n = 50
-- coordinates
x = linspace (n+1) (0,1::Double)
h = x @> 1
-- inner points
xInner = subVector 1 (n-1) x
mesh = [ (xj, yj) | xj <- toList xInner , yj <- toList xInner ]
-- laplace matrix
a = buildMatrix ((n-1)^2) ((n-1)^2) (idxf (n-1)) :: Matrix Double
-- right hand side
f = fromList $ map (\(x, y) -> (h^2)*x*(y^4)) mesh
-- solve linear system
loes = linearSolve a (asColumn f)
(me, loes, f) = poisson
```
%% Cell type:code id: tags:
``` haskell
import Graphics.Rendering.Plot
setPlots 1 1 >> withPlot (1, 1) (setDataset loes)
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` haskell
-- write out matrices to visualize with python
(xl,yl) = unzip me
xmat = reshape 49 (fromList xl)
ymat = reshape 49 (fromList yl)
-- create x and y matrix
saveMatrix "xmat.txt" "%f" xmat
saveMatrix "ymat.txt" "%f" ymat
saveMatrix "zmat.txt" "%f" loes
-- ../vis.py -x xmat.txt -y ymat.txt -z zmat.txt -t image
```
%% Cell type:markdown id: tags:
<img src="../images/poisson.png"></img>
%% Cell type:markdown id: tags:
## parallization
lets do the same in parallel using repa.
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## conjugate gradient solver
CG is a simple, efficient iterative solver for linear systems $A x = y$ with positive definite matrix A (i.e. $x^T A x \geq 0$ for all $x$).
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}
import Control.Monad.Identity
import Data.Array.Repa as R hiding ((++))
type Arr sh r = Array r sh Double
data CGState sh = CGState { cgx :: Arr sh U
, cgp :: Arr sh U
, cgr :: Arr sh U
, cgr2 :: Double
}
```
%% Cell type:markdown id: tags:
`Arr` is a type alias for convenience, `CGState` is the state of a CG iteration.
CG takes a function implementing a linear operator (type `Arr sh U -> Arr sh D`), a right hand side and an initial guess and returns the (lazy, infinte) list of the `CGState`s of all iterations.
Computations use Repa's parallel computation in a straight-forward way and needs to be performed in a monad. We use the `Identity` monad which does essentially nothing except for providing the sequencing. All operations are performed strictly using `BangPatterns`.
%% Cell type:code id: tags:
``` haskell
cg :: Shape sh => (Arr sh U -> Arr sh D) -> Arr sh U -> Arr sh U -> [CGState sh]
cg op rhs initial =
runIdentity $ do
!rInit <- computeP $ rhs -^ op initial
!r2Init <- normSquaredP rInit
return $ iterate cgStep (CGState initial rInit rInit r2Init)
where normSquaredP = sumAllP . R.map (^(2::Int))
scale a = R.map (* a)
-- the main iteration; maps a CGState to another CGState
cgStep (CGState x p r r2) =
runIdentity $ do
-- need inline type annotation here because computeP is polymorphic in its
-- output type (syntax like this requires ScopedTypeVariables)
!(q :: Arr sh U) <- computeP $ op p
!qp <- sumAllP $ q *^ p
let alpha = r2 / qp
!x' <- computeP $ x +^ scale alpha p
!r' <- computeP $ r -^ scale alpha q
!r2' <- normSquaredP r'
let beta = r2' / r2
!p' <- computeP $ r' +^ scale beta p
return $ CGState x' p' r' r2'
```
%% Cell type:markdown id: tags:
Using lazyness, the iterations can be decoupled from stopping rule and output of information.
%% Cell type:code id: tags:
``` haskell
takeUntil :: (a -> Bool) -> [a] -> [a]
takeUntil _ [] = []
takeUntil predicate (x:xs)
| predicate x = [x]
| otherwise = x : takeUntil predicate xs
process :: Monad m => [a] -> (a -> m ()) -> m a
process xs f = foldM (\_ x -> f x >> return x) undefined xs
runCG :: Shape sh => Double -> (Arr sh U -> Arr sh D) -> Arr sh U -> IO (Arr sh U)
runCG tol op rhs = do
let initial = computeS $ fromFunction (extent rhs) (const 0)
let steps' = cg op rhs initial
let r20 = cgr2 $ head steps'
let steps = takeUntil (\x -> sqrt (cgr2 x / r20) < tol) steps'
result <- process (zip [(1::Int)..] steps) $ \(n, cgs) ->
putStrLn $ show n ++ " " ++ show (sqrt $ cgr2 cgs / r20)
return $ cgx . snd $ result
```
%% Cell type:markdown id: tags:
- `takeUntil` is like `takeWhile`, but also returns the final iterate (why waste it?).
- `process` runs over a list, performs a monadic action on each element and returns the last element. It does not retain the start of the list while iterating, so it gets garbage collected properly.
- `runCG` executes `cg` with initial guess 0, a stopping rule based on residuals relative to the first initial guess, and outputs the residual while iterating.
### Example
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE FlexibleContexts #-}
import Data.Array.Repa.Stencil.Dim2
stencilOp :: Arr DIM2 U -> Arr DIM2 D
stencilOp = delay . mapStencil2 (BoundConst 0) st
where st = [stencil2| 0 1 0
1 0 1
0 1 0 |]
arr :: Arr DIM2 U
arr = computeS $ fromFunction (Z:.50:.50) $ \(Z:.i:.j) -> fromIntegral $ i+j
rhs :: Arr DIM2 U
rhs = computeS $ stencilOp arr
normS :: (Source r Double, Shape sh) => Arr sh r -> Double
normS = sqrt . sumAllS . R.map (^(2::Int))
do solution <- runCG 1e-3 stencilOp rhs
print $ normS (arr -^ solution) / normS arr
```
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment