Commit 7bc10581 authored by Jochen Schulz's avatar Jochen Schulz
Browse files

problem 2

parent 7f3d0359
# vim:ft=todo
LECTURE:
-exercise 23 : stiff nonstiff
-exercise 24: mandelbrot-irgendwas
-exercise 24: fft tiefpassfilter, jacobi mit repa, mandelbrot-irgendwas, poisson 1d
- jacobi poisson mit stencil
vorlesung 25: stencil muss auch erklaert werden (anhand differenzenquotient)
1d poisson beispiel fuer cg-loeser
25 poisson equation (2d statisch) mit repa (stencil!)
(-fft cg)
exercise 25: poisson 2d repa
26 continous optimization (minimierung: russel? siam -probleme )
projektionsverfahren (russell) sudoku-solver (bfgs fertig machen)
......@@ -31,6 +28,8 @@ TODO 2014-09-04 Musterloesungen erstellen
19 needs Musterloesungen
loesung problem 22.2
loesung 21.4
loesung 23.1 23.3
loesung 24.2
TODO 2015-08-14 change <img> to ihaskell-builtin types (because nb-viewer isnt working with this)
......@@ -54,6 +53,8 @@ xx *gerlind plonka waveleti bildverarbeitung
?? pde> meshless methods (heat equation, euler equation?
TODO 2015-09-01 data analysis testen (fehlen auf jeden fall pakete)
TODO 2015-08-05 allgemein: interpolation/fitting
......@@ -70,6 +71,7 @@ WAITING 2015-07-23 cloudhaskell
WAITING: 2015-07-23 11:08:19
mehr heuristische erklaerungen warum und wie etwas funktioniert.
---
TODO erklärung zu performance und polymorphie
......@@ -102,6 +104,8 @@ WAITING 2014-08-08 (Kategorientheorie?)
WAITING: 2015-06-12 12:28:13
ORGANISATION:
TODO 2014-09-04 evtl. eigenes Modul erstellen? oder noch andere mit implementieren?
TODO 2014-08-07 edgar onea kontaktieren https://sites.google.com/site/edgaroneagerman/
......
%% Cell type:markdown id: tags:
# Exercises 25
##problem 1
Put together the stationary heat equation (poisson equation) with the explicit Euler method to be able to solve the time dependent heat equation.
Given a rectangular domain $\Omega \subset \mathbb{R}^2 $ and a time dependent function
$u(x,t) , x \in \Omega , t \in \mathbb{R}^{+}$. Then the time dependent heat equation is
$$\frac{\partial u}{\partial t} - \alpha \triangle u = 0 \text{ in } \Omega $$
where $\alpha \in \mathbb{R}$ is a constant . Given dirichlet boundary conditions
$$u = R , \text{ on } \partial \Omega$$
with a function $R: \partial \Omega \mapsto C(\partial \Omega)$. At time $t=0$ let
$$u(x,0) = f(x), \forall x \in \Omega.$$
with a arbitrary, but fixed initial function $f: \mathbb{R} \mapsto \mathbb{R}$.
%% Cell type:code id: tags:
``` haskell
:extension NoMonomorphismRestriction
:ext BangPatterns
import Numeric.Container
import Numeric.LinearAlgebra.Algorithms
import Numeric.LinearAlgebra.Util
import Numeric.LinearAlgebra.HMatrix
import Graphics.Rendering.Plot
-- 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
al = 0.5
n = 50
poisson :: Int -> Double -> ([(Double, Double)],Matrix Double, Vector Double, Double)
poisson n al = (mesh, a, f, ht)
where
-- 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
-- For stability: calculate time-step:
ht = 2*h^2/( 4*al*h^2 )/100
(me,a,f,ht) = poisson n al
initialState :: Int -> Matrix Double
initialState n = buildMatrix (n-1) (n-1) fu
where
fu (i,j)
| (i - half)^2 + (j - half)^2 <= 4*radius^2 && (i - half)^2+(j - half)^2 >=radius^2 = 1
| otherwise = 0
where
half = div n 2
radius = div n 5
```
%% Cell type:code id: tags:
``` haskell
iter n f !x
| n==0 = x
| otherwise = iter (n-1) f (f x)
loes = iter 40 (\x -> x + (ht * al) `scale` (a #> x + f) ) (flatten $ initialState n)
--setPlots 1 1 >> withPlot (1, 1) (setDataset $ reshape (n-1) $loes )
```
%% Cell type:code id: tags:
``` haskell
-- write out matrices to visualize with python
(xl,yl) = unzip me
xmat = reshape (n-1) (fromList xl)
ymat = reshape (n-1) (fromList yl)
-- create x and y matrix
saveMatrix "xmat.txt" "%f" xmat
saveMatrix "ymat.txt" "%f" ymat
saveMatrix "zmat.txt" "%f" (reshape (n-1) loes)
-- ../vis.py -x xmat.txt -y ymat.txt -z zmat.txt -b "" -t image
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
## problem 2
Lets look into a bit image transformation. Take this greyscale image and make a *edge detection* on it. *edge detection* works through applying a convolution with
$$\left(\begin{array}
-- 1 & - 1 &- 1 \\
- 1 &8 &- 1 \\
- 1 &- 1 &- 1 \end{array}\right)$$
readImageFromBMP, mapStencil, repatomatrix, plot, graustufenbild angeben
%% Cell type:code id: tags:
``` haskell
import Data.Array.Repa.IO.BMP (readImageFromBMP, writeImageToBMP)
import Data.Array.Repa as R hiding ((++))
import Data.Array.Repa.Stencil
import Data.Array.Repa.Stencil.Dim2
import Data.Array.Repa.Algorithms.Pixel (doubleLuminanceOfRGB8, rgb8OfGreyDouble)
loadImage :: FilePath -> IO (Array U DIM2 Double)
loadImage file = do
loaded <- readImageFromBMP file
case loaded of
Left err -> error $ show err
Right img -> computeP $ R.map doubleLuminanceOfRGB8 img
saveImage :: FilePath -> Array U DIM2 Double -> IO ()
saveImage file img = do
!m <- foldAllP max 0 img
!img' <- computeP $ R.map (rgb8OfGreyDouble . (/ m)) img
writeImageToBMP file img'
if ( ( (x[i]-0.5)**2+(x[j]-0.5)**2 <= 0.1)
& ((x[i]-0.5)**2+(x[j]-0.5)**2>=.05) ):
u[i,j] = 1
{-# LANGUAGE QuasiQuotes #-}
edgeStencil :: Stencil DIM2 Double
edgeStencil = [stencil2| -1 -1 -1
-1 8 -1
-1 -1 -1 |]
img <- loadImage "../images/Grayscale_8bits_palette.bmp"
imgedge = computeS $ delay $ mapStencil2 (BoundConst 0) edgeStencil img
saveImage "../images/Grayscale_8bits_palette_edge.bmp" imgedge
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
## problem 3
- jacobi poisson mit stencil
%% Cell type:code id: tags:
``` haskell
-- For stability: calculate time-step:
ht = h**2*h**2/( 2*a*(h**2+h**2) )
-- Data for each time-step
data[t] + ht*a*(A*data[t] + b)
```
......
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