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

update lecture numeric part

parent fe8997f2
......@@ -14,48 +14,43 @@ TODO 2015-06-09 uebertragung der folien nach ihaskell notebooks und ins englisch
englisch fuer lecture 1
TODO 2015-06-26 install profiling libraries
cassava, htrace, monad-memo, monad-state
cassava, htrace, monad-memo, monad-state, hmatrix-gsl
TODO 2014-09-04 Musterloesungen erstellen
TODO 2015-06-24 cassava in 15 io: Data.csv ersetzen?
TODO 2015-06-24 anfang kompakter und kuerzer machen
TODO 2015-07-02 struktur vorlesung:
Mehr Beispiele Monaden2
Monaden2: State-Monade-Beispiele mit verschiedenen Typen
*hmatrix anfangen
*lenses (vor Chart!)
*charts
simple beispiele
ode poisson equation
lorentzattraktor
*ode: simple (lorentzattraktor)
ode: runge-kutta
*ode: boundary value problem poisson equation (1d)
integration
runge-kutta
repa parallelisierung poisson
interpolationi/fitting
ptimierung beispiele diskrete optimierung
optimierung beispiele diskrete optimierung
projektionsverfahren (russell) sudoku-solver
differentialgeometrie -> henrik
inverse probleme: intergralgleichungen,
gerlind plonka waveleti *bildverarbeitung
*gerlind plonka waveleti bildverarbeitung
pde> meshless methods (heat equation, euler equation?
pde zeitentwicklung *explicit euler
lube: finite elements (schwache formulierung(
diagrams evtl. weglassen
parallel and concurrent
cloudhaskell
TODO 2014-09-04 poisson: da ist noch ein Fehler drin
TODO 2014-09-02 sparse matrizen?
TODO 2014-08-21 immutable/mutable: was ist gemeint?
TODO 2015-06-12 monaden2 evtl. direkt hinter monaden1
TODO 2014-09-04 Mehr Beispiele Monaden2
TODO 2014-09-04 Monaden2: State-Monade-Beispiele mit verschiedenen Typen
---
......
......@@ -168,3 +168,25 @@ CANCELLED 2015-06-11 codemirror-standard languange setzen koennen
INHALT:
DONE 2015-06-24 cassava in 15 io: Data.csv ersetzen?
CLOSED: 2015-07-10 16:51:48
:LOGBOOK:
DONE: 2015-07-10 16:51:48
WAITING: 2015-07-10 16:51:48
DONE 2015-06-24 anfang kompakter und kuerzer machen
CLOSED: 2015-07-10 16:51:45
:LOGBOOK:
DONE: 2015-07-10 16:51:45
WAITING: 2015-07-10 16:51:44
DONE 2014-09-04 poisson: da ist noch ein Fehler drin
CLOSED: 2015-07-10 16:43:54
:LOGBOOK:
DONE: 2015-07-10 16:43:54
WAITING: 2015-07-10 16:43:53
DONE 2015-06-12 monaden2 evtl. direkt hinter monaden1
CLOSED: 2015-07-10 16:43:59
:LOGBOOK:
DONE: 2015-07-10 16:43:59
WAITING: 2015-07-10 16:43:59
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
%% Cell type:markdown id: tags:
#Poisson equation
# heat equation
%% Cell type:markdown id: tags:
## poisson equation in 2d
- 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)$.
**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
-- pure hmatrix
: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
--mwrap :: Matrix Double -> Int -> Int -> AlphaColour Double
--mwrap m x y = opaque ( rgb col col col )
-- where col = (m @@> (x,y))
-- do plotting with "plot" . used here but not explained further
fig = do
setPlots 1 1
withPlot (1,1) $ do
setDataset loes
addAxis XAxis (Side Lower) $ setTickLabelFormat DefaultTickFormat
addAxis YAxis (Side Lower) $ setTickLabelFormat DefaultTickFormat
setRangeFromData XAxis Lower Linear
setRangeFromData YAxis Lower Linear
fig
--b :: Diagram B
--b= rasterDia (mwrap loes) 700 700 -- # sized (dims2D 2 2)
```
%% Cell type:code id: tags:
``` haskell
import Graphics.Rendering.Plot
setPlots 1 1 >> withPlot (1, 1) (setDataset loes)
```
%% Output
......
%% Cell type:markdown id: tags:
# Numerik I
## Numerische Integration
## numeric integration
%% Cell type:code id: tags:
``` haskell
import Numeric.GSL.Integration
```
%% Cell type:markdown id: tags:
integrateQNG
: non-adaptive Gauss-Kronrod Integration im Intervall $(a,b)$.
**integrateQNG** <br /> non-adaptive Gauss-Kronrod integration in intervall $(a,b)$.
TODO: Erklaerung der Integration
TODO: Beispiele fuer eigens geschrieben Integration
%% Cell type:code id: tags:
``` haskell
integrateQNG :: Double -- Genauigkeit
-> (Double -> Double) --Funktion
-> Double -- a
-> Double -- b
-> (Double, Double) -- Ergebnis,Fehler
```haskell
integrateQNG
:: Double -- resolution
-> (Double -> Double) -- function
-> Double -- a
-> Double -- b
-> (Double, Double) -- result, error
```
%% Cell type:markdown id: tags:
*Beispiel*: $f(x) = \int_0^1 \frac{4}{1+x^2}$
*example*:
$f(x) = \int_0^1 \frac{4}{1+x^2}$
%% Cell type:code id: tags:
``` haskell
>integrateQNG 1E-6 (\x -> 4/(1+x^2)) 0 1
(3.1415926535897936,3.4878684980086326e-14)
integrateQNG 1E-6 (\x -> 4/(1+x^2)) 0 1
```
%% Output
%% Cell type:markdown id: tags:
## Numerische Differentation
## numeric differentation
%% Cell type:code id: tags:
``` haskell
import Numeric.GSL.Differentiation
```
%% Cell type:markdown id: tags:
derivCentral
: Differentation mit zentralen Differenzenquotient
**derivCentral** <br /> differentation with central difference quotient
%% Cell type:code id: tags:
``` haskell
derivCentral :: Double -- Schrittweite
-> (Double -> Double) -- Funktion
-> Double -- Ableitungspunkt
-> (Double, Double) -- Ergebnis und abs. Fehler
```haskell