Anged exercise with central difference quotient

TODO

@@ -20,12 +20,11 @@ DONE 2016-03-18 rewrite newton example with repa
 
 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) (pending..)
   ex24p3 solution missing
   ex24   mandelbrot-irgendwas 
   ex24 change to endtime?
+  ex25p1 zu repa konvertieren
   ex26p2 fehlt?
   ex27 inverse: weitere aufgaben?
   (-fft cg) aufgabe erstellen 

exercises/21 integration and differentiation.ipynb

 %% 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
 ```
 
 %% 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
 ```
 
 %% 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$.
+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 different central difference coefficents with a fixed $\sigma$.
+
+The simplest coefficents is
+$$\frac{\frac{-1}{2} f_{x-1} + \frac{1}{2} f_{x+1}}{h}$$
+
+which is expressed in coefficents (central is 0, and left and right the index goes down and up)
+
+$-1/2, 0,  1/2$
+
+Implement also the coefficients:
+
+$1/12, -2/3, 0,  2/3, -1/12$
+
+and
+
+$-1/60, 3/20, -3/4, 0, 3/4, -3/20, 1/60$
 
 *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
 ```
 
 %% 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
 ```

exercises/solutions/21 integration and differentiation.ipynb

 %% 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$.
+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 different central difference coefficents with a fixed $\sigma$.
+
+The simplest coefficents is
+$$\frac{\frac{-1}{2} f_{x-1} + \frac{1}{2} f_{x+1}}{h}$$
+
+which is expressed in coefficents (central is 0, and left and right the index goes down and up)
+
+$-1/2, 0,  1/2$
+
+Implement also the coefficients:
+
+$1/12, -2/3, 0,  2/3, -1/12$
+
+and
+
+$-1/60, 3/20, -3/4, 0, 3/4, -3/20, 1/60$
 
 *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
+periodicCentralDiff 1 (vec [1..12]) :: Vec Double
+```
+
+%% Cell type:code id: tags:
+
+``` haskell
+centralDiff :: Int -> Double -> Vec Double -> Vec Double
+centralDiff order h f = computeS $  fromFunction (Z:.n-6) (\(Z:.i) -> sumAllS $ scale (1/h) $ subVector i 7 f *^ diffp)
+    where
+        n = vlength f
+        diffp = case order of
+                1 -> vec [0, 0, -1/2, 0,  1/2, 0, 0] :: Vec Double
+                2 -> vec [0, 1/12, -2/3, 0,  2/3, -1/12, 0 ] :: Vec Double
+                3 -> vec [-1/60, 3/20, -3/4, 0, 3/4, -3/20, 1/60] :: Vec Double
 ```
 
 %% Cell type:code id: tags:
 
 ``` haskell
 normInf :: Vec Double -> Double
 normInf = maximum . map abs . toList
 
 -- 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))
+err :: Int -> Double -> Int -> Seed -> Double      -- one norm error for one seed
+err steps variance order seed = sqrt h * normInf (computeS $ subVector 3 (steps-3) 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))
+        fdn = centralDiff order h (computeS (ft +^ noise))
+        --fdn = periodicCentralDiff h (computeS (ft +^ noise))
 
 
-
-meanErr :: Int -> Int -> Double -> Double
-meanErr n steps variance = errs / fromIntegral n
+-- mean Error over several variances
+meanErr :: Int -> Int -> Double -> Int -> Double
+meanErr n steps variance order = errs / fromIntegral n
     where
-        errs = sum $ map (err steps variance) [0..n]
+        errs = sum $ map (err steps variance order) [0..n]
 
 
 
-makeResultMatrix :: [Int] -> [Double] -> Mat Double
-makeResultMatrix s v = computeS $ fromFunction (ix2 ns nv) $ run
+makeResultMatrix :: [Int] -> [Double] -> [Int] -> Mat Double
+makeResultMatrix s v o = computeS $ fromFunction (ix2 ns no) $ 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 50 (s!!i) (v!!j)
+          no = length o  -- number of orders
+          --run (Z :. i :. j) = meanErr 5 (s!!i) (v!!j)
+          run (Z :. i :. j) = meanErr 10 (s!!i) (v!!1) (o!!j)
 ```
 
 %% Cell type:code id: tags:
 
 ``` haskell
-res = makeResultMatrix [50,60..200] [0,0.2..1]
+-- import Graphics.Rendering.Plot
+res = makeResultMatrix [100,120..400] [0,0.0001..1] [1,2,3]
 res
--- toList res
-setPlots 1 1 >> withPlot (1, 1) (setDataset $ fromRepa res)
+-- setPlots 1 1 >> withPlot (1, 1) (setDataset $ fromRepa res)
+```
+
+%% Cell type:code id: tags:
+
+``` haskell
+import Graphics.Rendering.Chart.Easy
+import Graphics.Rendering.Chart.Backend.Cairo
+
+toRenderable $ Graphics.Rendering.Chart.Easy.plot (Graphics.Rendering.Chart.Easy.line "am" [zip [1..16] (toList $ column 1 res) ])
+toRenderable $ Graphics.Rendering.Chart.Easy.plot (Graphics.Rendering.Chart.Easy.line "am" [zip [1..16] (toList $ column 0 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.
+-- The plot shows -- as expected -- at first decreasing error with increasing number of steps. But too many steps
+-- make it worse: at some point the error *increases* with more points because making a too fine derivative is too
+-- dependant to the noise-level.
+-- --> There is a -- noise-dependent -- optimal step size.
 ```
 
 %% 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
 ```