Commit c562b1e5 by Max Voit

### practical feyerabend

parent b5fdb4d1
 %% Cell type:markdown id: tags: # boundary value problem (poisson equation in 1d) also known as a 2-point boundary value problem Find a function $$u:[0,1] \quad \rightarrow \quad \mathbb{R},$$ such that $$-u''(x) = e^x, \quad x \in (0,1)$$ $$u(0) = u(1) =0$$ There is no general closed solution because of the boundary conditions
$\Rightarrow$ numerical approximation of the solution. ## finite difference - discretice: $0=x_0 < \dots < x_{n}=1$ with $x_i=\frac{i}{n}$ - difference quotient: $$u''(x_i) \sim \frac{u(x_{i-1}) - 2 u(x_i) + u(x_{i+1})}{ h^2}, \quad h:=\frac{1}{n}$$ - plug in $-u''(x)=e^x$ yields $$-u(x_{i-1}) + 2 u(x_i) - u(x_{i+1}) = h^2 e^{x_i}, \quad i=1,\dots ,n-1$$ boundary condition $\Rightarrow$ $u(x_0)=u(x_n)=0$. These would lead to zero entries in the matrix and right hand side. - $\Rightarrow$ linear system for $u(x_1), \dots ,u(x_{n-1})$. Define $z=(z_1,\dots ,z_{n-1})^t=(u(x_1), \dots ,u(x_{n-1}))^t$. solve the linear system $Az=F$ with $$A:= \left( \begin{array} {ccccccc} 2 & -1 & & & 0 \\ -1 & 2 & -1 & & \\ & \ddots & \ddots & \ddots &\\ & & -1 & 2 & -1 \\ 0 & & & -1 & 2 \\ \end{array} \right), \ F:= h^2 \left( \begin{array}{c} e^\frac{1}{n}\\ \vdots \\ e^\frac{n-1}{n} \end{array} \right) .$$ ### Implementation - OLD %% Cell type:code id: tags:  haskell import Numeric.Container import Graphics.Rendering.Chart.Easy hiding (Matrix) import Graphics.Rendering.Chart.Backend.Cairo -- indexfunction for creating the matrix idxf (i,j) | i == j = 2 | i == j-1 || i == j+1 = -1 | otherwise = 0 twopoint = (x,a <\> (asColumn f)) -- solve the linear system where n = 100 -- divide the interval x = linspace n (0,1::Double) -- cut of boundary points (they are zero anyway) xi = subVector 1 (n-2) x -- build the matrix a = buildMatrix (n-2) (n-2) (idxf) :: Matrix Double -- right hand side f=exp(x) f = mapVector (\xi -> (1/(fromIntegral n))^2*exp(xi)) xi (x,sol) = twopoint -- fill up boundary values solfull = [0] ++ toList (flatten sol) ++ [0] toRenderable $do layout_title .= "2-point boundary value problem" plot (line "sol" [zip (toList x) solfull])  %% Cell type:markdown id: tags: ### HELPER %% Cell type:code id: tags:  haskell :ext FlexibleContexts import Data.Array.Repa import Numeric.LinearAlgebra.Helpers diag :: Num a => a -> Int -> Array D DIM2 a diag x n = fromFunction (ix2 n n) f where f (Z :. i :. j) | i == j = x | otherwise = 0 :t diag ident = diag 1 :t ident -- TODO check this computeS$ ident 3 :: Array U DIM2 Int  %% Output %% Cell type:code id: tags:  haskell :ext FlexibleContexts import Data.Array.Repa band :: Num a => [a] -> Int -> Array D DIM2 a band xs n = if even (length xs) then error "band: only odd vector lengths allowed!" else fromFunction (ix2 n n) f where m = length xs div 2 f (Z :. i :. j) | i == j = xs !! m | i < j && j-i <= m = xs !! (m+j-i) | j < i && i-j <= m = xs !! (m-i+j) | otherwise = 0 %% Cell type:code id: tags:  haskell computeS $band [1,2,3] 4 :: Array U DIM2 Int  %% Output %% Cell type:code id: tags:  haskell linspace :: (Fractional e) => Int -> (e, e) -> Array D DIM1 e linspace 0 _ = fromFunction (ix1 0) (const 0) linspace 1 (a,b) = fromFunction (ix1 1) (const (a+b/2)) linspace n (a,b) = fromFunction (ix1 n) f where f (Z :. x) = a + (b-a)*(fromIntegral x) / fromIntegral (n-1) (computeS$ linspace 5 (1,3)) :: Array U DIM1 Double  %% Output %% Cell type:markdown id: tags: ### Implementation %% Cell type:code id: tags:  haskell import Graphics.Rendering.Chart.Easy hiding (Matrix) import Graphics.Rendering.Chart.Backend.Cairo import Numeric.LinearAlgebra.Repa import Data.Array.Repa as R hiding ((++)) import Data.Array.Repa.Repr.ForeignPtr twopoint = (x,solveS a f) -- solve the linear system where n = 100 -- divide the interval x = linspace n (0,1::Double) -- cut of boundary points (they are zero anyway) xi = extract (ix1 1) (ix1 (n-2)) x -- build the matrix a = band [-1,2,-1] (n-2) -- right hand side f=exp(x) f = R.map (\xi -> (1/(fromIntegral n))^2*exp(xi)) xi (x,sol) = twopoint -- fill up boundary values solfull = [0] ++ toList sol ++ [0] toRenderable $do layout_title .= "2-point boundary value problem" plot (line "sol" [zip (toList x) solfull])  %% Output %% Cell type:code id: tags:  haskell  ... ... This diff is collapsed.  %% Cell type:markdown id: tags: # Runge-Kutta methods Runge-Kutta class methods are a implicit and explicit methods for temporal discretization for ODE's. Let it be a initial value problem $$y^\prime = f(t,y)$$ with$y(t_0)=y_0$. find$y(t)$. then runge-kutta methods have the form of $$y_{n + 1} = y_n + h \sum_{i = 1}^s b_i k_ i$$ $$k_i = f ( t_n + c_i h , y_n + h \sum_{j = 1}^s a_{ij} k_j )$$ with some initial value$y_0$and parameters$a_{ij}, b_i$and$c_i$## implicit Euler the simplest *implicit* runge-kutta method is the *implicit Euler method* (explicit is an exercise) with$c_i = 1$,$a_{ij} = 1$and$b_1 = 0$which then yields the update formula $$y_{n + 1} = y_n + h f ( t_{n+1} , y_{n+1} )$$ Now there is on both sides the unkown function and we have to solve an algebraic equation for this. This can be done with newton's method or with solving a linear equation depending on the problem. Because our example with the harmonic oscillator is linear but multidimensional we can solve a linear system of equations for$y_{n+1}\$ %% Cell type:code id: tags:  haskell import Numeric.Container import Data.Array.Repa