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

merge and work on exercises

parent d4d2f203
......@@ -22,7 +22,7 @@ TODO 2014-09-04 Musterloesungen erstellen
ex19 parallel arrays: rsa neu erstellen
ex20p5 jacobi solver
ex21p3 centralized differences + noise
ex22p4 auf cg umstellen (siam problem 20000)
ex22p4 auf cg umstellen (siam problem 20000) (pending..)
ex24p3 solution missing
ex24 mandelbrot-irgendwas
ex24 change to endtime?
......@@ -50,22 +50,25 @@ WAITING 2015-07-23 parallel and concurrent, cloudhaskell
:LOGBOOK:
WAITING: 2015-07-23 11:08:18
xx *gerlind plonka waveleti bildverarbeitung
?? pde> meshless methods (heat equation, euler equation?
xx lube: finite elements (schwache formulierung) heat equation (hmatrix)
xx+1 *differentialgeometrie -> henrik (laplace weak implementieren) (hmatrix)
xx *inverse probleme: intergralgleichungen,
xx *gerlind plonka waveleti bildverarbeitung
volterra intergralgleichungen (numeric differentiation) (schon in ex27?)
?? 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
TODO 2014-08-21 immutable/mutable
WAITING 2014-08-08 algebraische typen detaillierter
DONE 2014-08-08 algebraische typen detaillierter
CLOSED: 2016-03-21 16:13:21
:LOGBOOK:
DONE: 2016-03-21 16:13:21
WAITING: 2015-07-02 17:28:15
WAITING 2014-09-01 (spline-interpolation? , (siehe auch mapy))
......@@ -80,6 +83,8 @@ WAITING 2014-07-29 headmay : http://hackage.haskell.org/package/safe-0.3/docs/Sa
:LOGBOOK:
WAITING: 2015-06-12 12:28:11
ANDEREVORLESUNG:
WAITING 2014-08-05 (Lambda Kalkuel einfuehren)
:LOGBOOK:
WAITING: 2015-06-12 12:28:12
......
-- ghc -O2 -prof -fprof-auto -rtsopts factorial.hs && ./factorial +RTS -pi -hc -sstderr -RTS
{-# LANGUAGE BangPatterns #-}
facE :: (Show a, Integral a) => a -> a
facE i
| i<0 = error "Negative Faculty not defined"
......@@ -20,6 +21,8 @@ facRf x = fac' x 1 where
fac' 1 y = y
fac' f y = fac' (f-1) $! (f*y)
Fractional
main :: IO ()
main = do
let g = facNE 10000
......
%% Cell type:markdown id: tags:
# Exercises: type classes
## problem 1
Explain the difference between
```haskell
type Point = (Double, Double)
```
and
```haskell
data Point = Point Double Double
```
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 2
Write an instance of `Show` for `Shape`
```haskell
data Point = Point Double Double
instance Show Point where
show (Point x y) = "(" ++ show x ++ ", " ++ show y ++ ")"
data Shape = Circle Point Double
| Rectangle Point Double Double
```
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 3
- Write an instance of `Eq` for `Shape`.
- Write an instance of `Ord` for `Shape`.
- Compare to the behaviour of the automatically generated `Ord` instance from a `deriving` clause and explain.
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 4
Vectors (in the mathematical sense) are required in many numerical applications.
- Implement a typeclass for vectors of arbitrary dimensions.
- Implement addition, scalar multiplication and inner product.
- How could matrices be represented?
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 5
Write classes `ShowHex` and `ShowBin` which output numbers in hex and binary, and implement them for `Int`.
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 6
lookup the function `unfoldr` and first write an own version of `iterate` by using `unfoldr`. Then rewrite the instances of problem 5 via using `unfoldr`.
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 7
Converting between metric and imperial units: Write a function `convertUnit :: (Double, String) -> (Double, String)` that behaves like
```haskell
convertUnit (1, "m") == (1.09361, "yd")
convertUnit (1, "kg") == (2.20462, "lb")
```
etc...
Is there a more type-safe way of doing this?
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 8
Prove the second functor law for lists:
```haskell
fmap (f . g) == fmap f . fmap g
```
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 9
Create a data type for polynominial rings with real coefficients and write the instances for `Show`, `Eq` and `Num`.
Create a data type for polynomial rings with real coefficients and write the instances for `Show`, `Eq` and `Num`.
*Note*: a formal polynominial is nothing else than a list of coefficents.
*Note*: a formal polynomial is nothing else than a list of coefficents.
%% Cell type:code id: tags:
``` haskell
```
......
%% Cell type:markdown id: tags:
# Exercises: boundary value problems
## problem 1
implement a *Gauss-Seidel* algorithm for solving linear equations.
The Gauss–Seidel method is an iterative technique for solving a square system of n linear equations with unknown x:
$$A x = b$$
It is defined by the iteration
$$L_∗ x^{( k + 1 )} = b − U x^{( k )}$$
where $x^{( k )}$ is the $k$-th approximation or iteration of $x$ ,$x^{k + 1}$ is the next value of $x$ , and the matrix $A$ is decomposed into a lower triangular component $L_∗$ , and a strictly upper triangular component $U$: $A = L_∗ + U$
We rewrite this a little bit to for the algorithm:
$$x^{( k + 1 )} = T x^{( k )} + C$$
where
$$T = - L_∗^{-1} U$$ and $$C = L_∗^{-1} b$$
Write a recursive function which makes this iteration. The iteration should stop both from a tolerance value of the difference of adjacent $x$ and a maximum of iterations.
For a test case use the matrix and right hand side from the lecture.
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 2
Implement the *jacobi* iterative algorithm for solving linear equations. First decompose the
matrix $A$ into the diagonal part $D$ and the remnant $R$:
$$ A = D + R$$
Then iterate the following
$$x^{( k + 1 )} = D^{ - 1} ( b - R x^{( k )} )$$
where $x^{ ( k )}$ is the $k$-th approximation or iteration of $x$ and $x^{ ( k + 1 )}$ is the next.
Write a recursive function which makes this iteration. The iteration should stop both from a tolerance value of the difference of adjacent x and a maximum of iterations. For a testcase use the matrix and right hand side from the lecture.
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 3
Write an implementation of the *conjugate gradient* (CG) method. This method solves
$$A x = b$$
by an iteration of the form $x_0 = 0$, $r_0 = p_0 = b$,
$$
\begin{align}
\alpha_k &= \frac{{\lVert r_k \rVert}^2}{\langle p_k, A p_k \rangle} \\
x_{k+1} &= x_k + \alpha_k p_k \\
r_{k+1} &= r_k - \alpha_k A p_k \\
p_{k+1} &= r_{k+1} + \frac{{\lVert r_{k+1} \rVert}^2}{{\lVert r_k \rVert}^2} p_k
\end{align},
$$
stopping the iterations at the smallest $K$ for which $\lVert r_K \rVert < \text{tol}$.
An interesting property of this method in contrast to the ones from the previous problems is that it does not need $A$ in form of a matrix; it suffices to be able to apply $A$ to any given vector $p$, so $A$ can also be passed as a function `Vec Double -> Vec Double`.
The CG implementation should therefore have a signature like this:
```haskell
cg :: Double -> (Vec Double -> Vec Double) -> Vec Double -> Vec Double
cg tol f rhs = [...]
```
Apply the method to the same test case as above, but this time reimplement it using the function
```haskell
conv :: Vec t -> Vec t -> Vec t
```
from `Numeric.LinearAlgebra.Repa`, which implements a convolution without actually creating the matrix. Note that `conv` returns a vector that is larger than the input; use the boundary values to handle this.
%% Cell type:code id: tags:
``` haskell
```
%% Cell type:markdown id: tags:
## problem 4
Let A be a $20.000 \times 20.000$ matrix, which entries are all zero besides the prime numbers $2, 3, 5, 7, \ldots, 224737$ on the diagonal and the number 1 in all entries $a_{ij}$ with
$$|i − j| = 1, 2, 4, 8, \ldots , 16384.$$
What is the $(1, 1)$ entry of $A^{−1}$? How can you check/visualize the structure of the matrix? Be sure to use a good method for saving the matrix (use sparse..).
*bonus*: implement sparse variants of gauss-seidel and/or conjugate gradiant and try to use them as solver here.
What is the $(1, 1)$ entry of $A^{−1}$? How can you check/visualize the structure of the matrix? Be sure to use a good method for saving the matrix.
%% Cell type:code id: tags:
``` haskell
```
......
%% Cell type:markdown id: tags:
# Exercises: type classes
## problem 1
Explain the difference between
```haskell
type Point = (Double, Double)
```
and
```haskell
data Point = Point Double Double
```
%% Cell type:markdown id: tags:
## problem 2
Write an instance of `Show` for `Shape`
```haskell
data Point = Point Double Double
instance Show Point where
show (Point x y) = "(" ++ show x ++ ", " ++ show y ++ ")"
data Shape = Circle Point Double
| Rectangle Point Double Double
```
%% Cell type:code id: tags:
``` haskell
data Point = Point Double Double
instance Show Point where
show (Point x y) = "(" ++ show x ++ ", " ++ show y ++ ")"
data Shape = Circle Point Double
| Rectangle Point Double Double
instance Show Shape where
show (Circle c r) = "Circle with center " ++ show c ++ " and radius " ++ show r
show (Rectangle p a b) = "Rectangle with top-left corner " ++ show p ++
" and edges " ++ show a ++ " and " ++ show b
```
%% Cell type:markdown id: tags:
## problem 3
- Write an instance of `Eq` for `Shape`.
- Write an instance of `Ord` for `Shape`.
- Compare to the behaviour of the automatically generated `Ord` instance from a `deriving` clause and explain.
%% Cell type:code id: tags:
``` haskell
instance Eq Point where
(Point x1 y1) == (Point x2 y2) = (x1 == x2) && (y1 == y2)
instance Eq Shape where
(Circle c1 r1) == (Circle c2 r2) = (c1 == c2) && (r1 == r2)
(Rectangle p1 a1 b1) == (Rectangle p2 a2 b2) = (p1 == p2) && (a1 == a2) && (b1 == b2)
_ == _ = False
```
%% Cell type:code id: tags:
``` haskell
import Data.Function
-- problem: unclear how to compare shapes. e.g. compare areas:
instance Ord Shape where
compare = compare `on` area
where area (Rectangle _ a b) = a * b
area (Circle _ r) = pi * r**2
```
%% Cell type:code id: tags:
``` haskell
data PointDeriv = PointDeriv Double Double deriving (Eq, Ord)
data ShapeDeriv = CircleDeriv PointDeriv Double
| RectangleDeriv PointDeriv Double Double
deriving (Eq, Ord)
```
%% Cell type:code id: tags:
``` haskell
CircleDeriv (PointDeriv 0 0) 1 < CircleDeriv (PointDeriv 1 0) 10
```
%% Cell type:code id: tags:
``` haskell
CircleDeriv (PointDeriv 0 0) 10 < RectangleDeriv (PointDeriv 0 0) 1 1
```
%% Cell type:code id: tags:
``` haskell
RectangleDeriv (PointDeriv 0 0) 2 3 < RectangleDeriv (PointDeriv 0 0) 3 2
```
%% Cell type:markdown id: tags:
## problem 4
Vectors (in the mathematical sense) are required in many numerical applications.
- Implement a typeclass for vectors of arbitrary dimensions.
- Implement addition, scalar multiplication and inner product.
- How could matrices be represented?
%% Cell type:code id: tags:
``` haskell
data Vector = Vector { vec :: [Double] } deriving (Show)
-- Downside: length of vector is unspecified, needs to be checked at runtime.
-- Otherwise, we would need a seperate type for each dimension: Vector1d, Vector2d, ...
addVector :: Vector -> Vector -> Vector
addVector a b
| length c /= length d = error "Dimensions are not equal"
| otherwise = Vector (zipWith (+) c d)
where c = vec a
d = vec b
scalar :: Vector -> Double -> Vector
scalar a b = Vector (map (*b) (vec a))
scalarProd :: Vector -> Vector -> Double
scalarProd a b
| length c /= length d = error "Vector sizes must match"
| otherwise = sum (zipWith (*) c d)
where c = vec a
d = vec b
```
%% Cell type:code id: tags:
``` haskell
data Matrix = Matrix { columns :: [Vector] }
```
%% Cell type:markdown id: tags:
## problem 5
Write classes `ShowHex` and `ShowBin` which output numbers in hex and binary, and implement them for `Int`.
%% Cell type:code id: tags:
``` haskell
class ShowHex a where
showHex :: a -> String
instance ShowHex Int where
showHex 0 = "0x0"
showHex n = "0x" ++ (reverse $ showHex' n)
where showHex' 0 = ""
showHex' m = let (q, r) = quotRem m 16
in ((['0'..'9']++['a'..'f']) !! r) : showHex' q
```
%% Cell type:code id: tags:
``` haskell
class ShowBin a where
showBin :: a -> String
instance ShowBin Int where
showBin 0 = "0"
showBin n = reverse $ showBin' n
where showBin' 0 = ""
showBin' m = let (q, r) = quotRem m 2
in (['0','1'] !! r) : showBin' q
```
%% Cell type:markdown id: tags:
## problem 6
lookup the function `unfoldr` and first write an own version of `iterate` by using `unfoldr`. Then rewrite the instances of problem 5 via using `unfoldr`.
%% Cell type:code id: tags:
``` haskell
import Data.List
iterate' f = unfoldr (\x -> Just (x, f x))
```
%% Cell type:code id: tags:
``` haskell
-- The recursion pattern above is abstracted by
-- unfoldr :: (a -> Maybe (b, a)) -> a -> [b]
-- So a nicer way to write the instances is this:
import Data.List
import Data.Tuple
toRadix rad 0 = [0]
toRadix rad n = reverse . unfoldr (\m -> if m == 0 then Nothing else Just (swap $ quotRem m rad)) $ n
toDigit = ((['0'..'9'] ++ ['a'..'z']) !!)
instance ShowHex Int where
showHex = ("0x" ++) . map toDigit . toRadix 16
instance ShowBin Int where
showBin = ("0x" ++) . map toDigit . toRadix 2
```
%% Cell type:markdown id: tags:
## problem 7
Converting between metric and imperial units: Write a function `convertUnit :: (Double, String) -> (Double, String)` that behaves like
```haskell
convertUnit (1, "m") == (1.09361, "yd")
convertUnit (1, "kg") == (2.20462, "lb")
```
etc...
Is there a more type-safe way of doing this?
%% Cell type:code id: tags:
``` haskell
convertUnit :: (Double, String) -> (Double, String)
convertUnit (x, l)
| l == "m" = (1.09361 * x, "yd")
| l == "l" = (0.264172 * x, "gal")
| l == "kg" = (2.20462 * x, "lb")
| l == "yd" = (x / 1.09361, "m")
| l == "gal" = (x / 0.264172, "l")
| l == "lb" = (x / 2.20462, "kg")
| otherwise = error "Invalid input"
```
%% Cell type:code id: tags:
``` haskell
-- Using types (via http://shuklan.com/haskell/lec05.html):
data MetricUnit = Meter
| Liter
| KiloGram
deriving (Show, Eq)
```
%% Cell type:code id: tags:
``` haskell
data ImperialUnit = Yard
| Gallon
| Pound
deriving (Show, Eq)
```
%% Cell type:code id: tags:
``` haskell
data Measurement = MetricMeasurement Double MetricUnit
| ImperialMeasurement Double ImperialUnit
deriving (Show, Eq)
```
%% Cell type:code id: tags:
``` haskell
convertUnit :: Measurement -> Measurement
convertUnit (MetricMeasurement x u)
| u == Meter = ImperialMeasurement (1.0936*x) Yard
| u == Liter = ImperialMeasurement (0.2642*x) Gallon
| u == KiloGram = ImperialMeasurement (2.2046*x) Pound
convertUnit (ImperialMeasurement x u)
| u == Yard = MetricMeasurement (0.9144*x) Meter
| u == Gallon = MetricMeasurement (3.7854*x) Liter
| u == Pound = MetricMeasurement (0.4536*x) KiloGram
```
%% Cell type:markdown id: tags:
## problem 8
Prove the second functor law for lists:
```haskell
fmap (f . g) == fmap f . fmap g
```
%% Cell type:code id: tags:
``` haskell
-- Induction start:
fmap (f . g) [] == [] == fmap f (fmap g [])
-- Induction step:
fmap (f . g) (x:xs)
== (f . g) x : fmap (f . g) xs
== f (g x) : fmap f (fmap g xs)
== fmap f (g x : fmap g xs)
== fmap f (fmap g (x:xs))
```
%% Cell type:markdown id: tags:
## problem 9
Create a data type for polynominial rings with real coefficients and write the instances for `Show`, `Eq` and `Num`.
Create a data type for polynomial rings with real coefficients and write the instances for `Show`, `Eq` and `Num`.
*Note*: a formal polynominial is nothing else than a list of coefficents.
*Note*: a formal polynomial is nothing else than a list of coefficents.
%% Cell type:code id: tags:
``` haskell
-- TODO : cleanup
import Data.List
newtype Polynomial = Polynomial [Double]
cleanPoly (Polynomial a) = Polynomial (dropWhileEnd (==0) a)
deg :: Polynomial -> Int
deg (Polynomial a) = length a - 1
lcoeff :: Polynomial -> Double
lcoeff (Polynomial a) = last a
monomial :: Int -> Double -> Polynomial
monomial d a = Polynomial (replicate d 0 ++ [a])
fromDouble :: Double -> Polynomial
fromDouble a = Polynomial [a] -- also: fromDouble = monomial 0
zero = Polynomial []
printPoly a = ((++) "P(X) = " . intercalate " + " . reverse . filter (/= "") . zipWith zipper [0..length a-1]) a
where
zipper _ 0 = ""
zipper power coeff = show coeff ++ "*X^" ++ show power
instance Show Polynomial where
show (Polynomial a) = printPoly a
instance Eq Polynomial where
(Polynomial a) == (Polynomial b) = a == b
addPoly [] p = p
addPoly p [] = p
addPoly (x:xs) (y:ys) = (x+y:addPoly xs ys)
cauchyProd p q = [sum [(p !! i) * (q !! j) | i <- [0..lp-1], j <- [0..lq-1], i+j == n] | n <- [0..lp+lq-2]]
where
lp = length p
lq = length q
instance Num Polynomial where
(Polynomial a) + (Polynomial b) = cleanPoly (Polynomial (addPoly a b))
negate (Polynomial a) = Polynomial (map negate a)
fromInteger a = Polynomial [fromInteger a]
abs (Polynomial a) = Polynomial (map abs a)
signum (Polynomial a) = Polynomial (map signum a)
(Polynomial a) * (Polynomial b) = Polynomial (cauchyProd a b)
```
%% Cell type:code id: tags:
``` haskell
-- example
p1 = fromDouble 3.14
p2 = monomial 3 2.17
p1
p2
```
%% Cell type:code id: tags:
``` haskell
```
......
%% Cell type:markdown id: tags:
# Exercises: integration and differentiation
## problem 1
Find an approximate solution for the fixed point equation
$$x = e^{-x}.$$
Use $x_0 = 1.1$ as initial value and break after 1000 iterations or if the error is small enough.
%% Cell type:code id: tags:
``` haskell
fixp :: Double -> Int -> Double -> (Double -> Double) -> Double
fixp tol maxit x0 f =
case converged of
x:_ -> x
otherwise -> last iterates
where iterates = take maxit $ iterate f x0
converged = [ x | (x, xprev) <- zip (tail iterates) iterates
, abs (x - xprev) < tol ]
fixp 1e-5 1000 1.1 $ \x -> exp (-x)
```
%% Cell type:markdown id: tags:
## problem 2
Write a function that performs simple numerical (Riemann) integration using step functions.
A function $\phi:[a,b] \rightarrow \mathbb{R}$ is called step function if there exists a partion of $[a,b]$ into intervals such that $\phi$ is constant on each interval.
%% Cell type:code id: tags:
``` haskell
int :: Double -> Double -> Int -> (Double -> Double) -> Double
int a b n f = (b-a) / n' * sum (map f xs)
where xs = [ a + fromIntegral k / n' * (b-a) | k <- [0..n] ]
n' = fromIntegral n
```
%% Cell type:markdown id: tags:
## problem 3
Numerically compute the derivative of $f(x) = \sin(15x)$, $x \in [0, 2\pi]$ given a vector of noisy values $f(x_k) + \epsilon_k$, for $x_k = 2 \pi k / N$, $k \in \{0, 1, \ldots, N-1\}$. Here, $\epsilon_k$ are normally distributed noise variables with mean 0 and variance $\sigma$. Compare the result to the analytical derivative by computing a suitably scaled norm of the difference, averaged over sufficiently many noise realizations, and present the results graphically for varying $N$ and $\sigma$.
*Note*: You can generate random vectors either using the standard `System.Random` functions, or
```haskell
randomVector :: Seed -> RandDist -> Int -> Vec Double
```
from `Numeric.LinearAlgebra.Repa`. As seed values for the different noise realization, you can use `[0..]` for simplicity.
%% Cell type:code id: tags:
``` haskell
import Graphics.Rendering.Plot
import qualified Data.Array.Repa as Repa
import Data.Array.Repa hiding (map, (++))
import Numeric.LinearAlgebra.Helpers
import Numeric.LinearAlgebra.Repa
```
%% Cell type:code id: tags:
``` haskell
periodicCentralDiff :: Double -> Vec Double -> Vec Double
periodicCentralDiff h f = computeS $ scale (1/(2*h)) (computeS (fLeft -^ fRight) :: Vec Double)
where
n = vlength f
fLeft = subVector 1 (n-1) f `append` subVector 0 1 f
fRight = subVector (n-1) 1 f `append` subVector 0 (n-1) f
periodicCentralDiff 1 (vec [1,2,3,2,1,0]) :: Vec Double
```
%% Cell type:code id: tags:
``` haskell
steps = 200
xs = subVector 0 steps $ computeS $ linspace (steps + 1) (0, 2*pi)
--xs
fs :: Vec Double
fs = computeS $ Repa.map (\x -> sin (15*x)) xs
--fs
dfs :: Vec Double
dfs = computeS $ Repa.map (\x -> 15*cos (15*x)) xs
--dfs
```
%% Cell type:code id: tags:
``` haskell
h = 2*pi/(steps+1)
--h
--periodicCentralDiff h fs
```
%% Cell type:code id: tags:
``` haskell
computeS $ Repa.map (\x -> sin (x*2)) (vec [0,1,3.14,6.28]) :: Vec Double
```
%% Cell type:code id: tags:
``` haskell
normInf :: Vec Double -> Double
normInf = maximum . map abs . toList
errs :: Int -> Double -> [Double]
errs steps variance = map (err . makeNoise) [0..]
where
xs :: Vec Double -- all x locations
xs = subVector 0 steps $ computeS $ linspace (steps + 1) (0, 2*pi)
f :: Vec Double -- f at all x locations
f = computeS $ Repa.map (\x -> sin (15*x)) xs
f'exact :: Vec Double -- exact derivative at the locations
f'exact = computeS $ Repa.map (\x -> 15*cos (15*x)) xs
h :: Double -- the step size
h = xs!(ix1 1) - xs!(ix1 0)
makeNoise :: Seed -> Vec Double -- generate noise for a given seed
makeNoise seed = sqrt variance `hscale` randomVector seed Gaussian steps
err :: Vec Double -> Double -- the accumulated error in all noisy derivatives
err noise = let f' = periodicCentralDiff h (computeS (f +^ noise))
in sqrt h * normInf (computeS (f'exact -^ f'))
-- TODO: this is quite haskell^Whaesslich
```
-- create vals x and step size h
coords :: Int -> (Vec Double, Double)
coords steps = (xs,h)
where
xs = subVector 0 steps $ computeS $ linspace (steps + 1) (0, 2*pi)
h = xs!(ix1 1) - xs!(ix1 0)
err :: Int -> Double -> Seed -> Double -- one norm error for one seed
err steps variance seed = sqrt h * normInf (computeS (fdt -^ fdn))
where
(xs,h) = coords steps
ft = computeS $ Repa.map (\x -> sin (15*x)) xs :: Vec Double
fdt = computeS $ Repa.map (\x -> 15*cos (15*x)) xs ::Vec Double-- true central difference
-- generate noise for a given seed
noise = sqrt variance `hscale` randomVector seed Gaussian steps
fdn = periodicCentralDiff h (computeS (ft +^ noise))
%% Cell type:code id: tags:
``` haskell
meanErr' :: Int -> Int -> Double -> Double
meanErr' n steps variance = (sum . take n $ errs steps variance) / fromIntegral n
meanErr :: Int -> Double -> Double
meanErr = meanErr' 100
meanErr :: Int -> Int -> Double -> Double
meanErr n steps variance = errs / fromIntegral n
where
errs = sum $ map (err steps variance) [0..n]
makeResultMatrix :: [Int] -> [Double] -> Mat Double
makeResultMatrix s v = computeS $ fromFunction (ix2 ns nv) $ run
where ns = length s -- number of steps
nv = length v -- number of different variances
-- we plot the square of the mean error for better contrast
run (Z :. i :. j) = meanErr (s!!i) (v!!j) ^ 2
run (Z :. i :. j) = meanErr 50 (s!!i) (v!!j)
```
res = makeResultMatrix [50,60..200] [0,0.2..5]
--toList res
%% Cell type:code id: tags:
``` haskell
res = makeResultMatrix [50,60..200] [0,0.2..1]
res
-- toList res
setPlots 1 1 >> withPlot (1, 1) (setDataset $ fromRepa res)
```
%% Cell type:code id: tags:
``` haskell
-- The plot shows -- as expected -- an increasing error with increasing noise
-- variance, but also a better robustness with larger step sizes. On the other
-- hand, too large step sizes lead to large errors since the approximation
-- becomes bad. There is a -- noise-dependent -- optimal step size somewhere
-- in the middle which balances these errors.
```
%% Cell type:markdown id: tags:
## problem 4
Consider the one-dimensional integral equation
$$\int_0^1 \phi(x)e^{-|x-y|}\sin(|x-y|)^2 dx = f(y),\quad y \in [0, 1]$$
where $\phi(x) \in \mathbb{R}$ and $f(y) = \sin(y)$. Approximate this by a discretized version
$$A \phi = f$$
with a matrix $A$ and solve the equation numerically. Check if the solution actually solves the discrete system.
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE FlexibleContexts #-}
import Data.Array.Repa
import Numeric.LinearAlgebra.Repa
import Numeric.LinearAlgebra.Helpers
kernel :: Double -> Double
kernel x = exp (-z) * sin z ^ (2::Int)
where z = abs x
(//) :: Integral a => a -> a -> Double
x // y = fromIntegral x / fromIntegral y
infixl 7 //
mkMatrix :: Int -> Int -> Mat Double
mkMatrix ny nx = mat
where mat = computeS $ fromFunction (ix2 ny nx) get
get (Z :. iy :. ix) = (1//nx) * kernel ((ix // nx) - (iy // ny))
rhs :: Int -> Vec Double
rhs n = computeS $ fromFunction (ix1 n) $ \(Z:.k) -> sin (k // n)
run nx ny = norm2 $ computeS (y -^ y')
where mat = mkMatrix ny nx
y = rhs ny
x = mat <\> y
y' = mat `app` x
run 100 100
```
%% Cell type:code id: tags:
``` haskell
```
......
%% Cell type:markdown id: tags:
# Exercises: boundary value problems
## problem 1
implement a *Gauss-Seidel* algorithm for solving linear equations.
The Gauss–Seidel method is an iterative technique for solving a square system of n linear equations with unknown x:
$$A x = b$$
It is defined by the iteration
$$L_∗ x^{( k + 1 )} = b − U x^{( k )}$$
where $x^{( k )}$ is the $k$-th approximation or iteration of $x$ ,$x^{k + 1}$ is the next value of $x$ , and the matrix $A$ is decomposed into a lower triangular component $L_∗$ , and a strictly upper triangular component $U$: $A = L_∗ + U$
We rewrite this a little bit to for the algorithm:
$$x^{( k + 1 )} = T x^{( k )} + C$$
where
$$T = - L_∗^{-1} U$$ and $$C = L_∗^{-1} b$$
Write a recursive function which makes this iteration. The iteration should stop both from a tolerance value of the difference of adjacent $x$ and a maximum of iterations.
For a test case use the matrix and right hand side from the lecture.
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE FlexibleContexts #-}
import Numeric.LinearAlgebra.Repa
import Numeric.LinearAlgebra.Helpers
import qualified Data.Array.Repa as Repa
import Data.Array.Repa hiding (map, (++))
-- lets first create a matrix a and right hand side as in the lecture
idxf (Z :. i :. j)
| i == j = 2
| i == j-1 || i == j+1 = -1
| otherwise = 0
setupsystem :: Int -> (Mat Double, Vec Double, Vec Double)
setupsystem n = (a,f,xi)
where
-- divide the interval
x :: Vec Double
x = computeS $ linspace n (0,1::Double)
-- cut of boundary points (they are zero anyway)
xi :: Vec Double
xi = subVector 1 98 x
-- build the matrix
a :: Mat Double
a = computeS $ fromFunction (ix2 (n-2) (n-2)) (idxf) :: Mat Double
f :: Vec Double
f = computeS $ Repa.map (\xi -> (1/(fromIntegral n))^2*exp(xi)) xi
```
%% Cell type:code id: tags:
``` haskell
gaussSeidel :: Mat Double -> Vec Double -> Int -> Double -> (Vec Double, Double)
gaussSeidel a b n tol = gsit (computeS $