src/FEM/R1.hs

   1 {-# LANGUAGE GADTs #-}
   2 {-# LANGUAGE ScopedTypeVariables #-}
   3
   4 -- | Finite Element Method (FEM) to solve the problem,
   5 --
   6 --     -(A(x)*u'(x))' + c(x)*u*(x) = f(x) on (a,b)
   7 --
   8 --   with one of the two types of boundary conditions,
   9 --
  10 --     * Dirichlet: u(a) = u(b) = 0
  11 --
  12 --     * A(a)u'(a) = alpha, A(b)u'(b) = beta
  13 --
  14 --   if c(x) = 0, then it is assumed that the boundary conditions are
  15 --   Dirichlet.
  16 --
  17 --   The code creates a linear system,
  18 --
  19 --     (K + M)x = F
  20 --
  21 --   to be solved for @x@, where @K@ (\"big_K\") is the stiffness
  22 --   matrix, @M@ (\"big_M\") is the mass matrix, and @F@ (\"big_F\")
  23 --   is the load vector.
  24 --
  25 --   Warning: until PLU factorization is implemented, we can only
  26 --   solve the resulting system if it's positive definite!
  27 --
  28 module FEM.R1
  29 where
  30
  31 import qualified Algebra.Algebraic as Algebraic ( C )
  32 import qualified Algebra.Field as Field ( C )
  33 import qualified Algebra.RealField as RealField ( C )
  34 import qualified Algebra.ToRational as ToRational ( C )
  35 import Data.Vector.Fixed ( Arity, S )
  36 import NumericPrelude
  37 import qualified Prelude as P
  38
  39 import Integration.Gaussian ( gaussian )
  40 import Linear.Matrix (
  41   Col,
  42   Mat(..),
  43   Row,
  44   (!!!),
  45   construct,
  46   element_sum2,
  47   fromList,
  48   ifoldl2,
  49   map2,
  50   nrows,
  51   rows2,
  52   set_idx,
  53   toList,
  54   transpose,
  55   zip2,
  56   zipwith2 )
  57 import Linear.System ( solve_positive_definite )
  58 import Piecewise ( Piecewise(..), from_intervals )
  59 import Polynomials.Orthogonal ( legendre )
  60
  61 -- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed,
  62 --   there's no additional information conveyed by this type.
  63 data Dirichlet a = Dirichlet { domain_dirichlet :: Interval a }
  64
  65 -- | Neumann boundary conditions. @alpha@ specifies A(a)u'(b) and
  66 --   @beta@ specifies A(b)u'(b).
  67 data Neumann a =
  68   Neumann { domain_neumann :: Interval a,
  69             alpha :: a,
  70             beta :: a }
  71
  72 -- | Boundary conditions can be either Dirichlet or Neumann.
  73 type BoundaryConditions a = Either (Dirichlet a) (Neumann a)
  74
  75 type Interval a = (a,a)
  76
  77 -- | The data for the PDE that we are attempting to solve.
  78 data PDE a =
  79   PDE {
  80     -- | A(x)
  81     big_A   :: (a -> a),
  82     -- | c(x)
  83     c      :: (a -> a),
  84     -- | f(x)
  85     f      :: (a -> a),
  86
  87     -- | The boundary conditions. The domain also specifies the
  88     --   boundary in R^1.
  89     bdy    :: BoundaryConditions a }
  90
  91
  92
  93 -- | Non-PDE parameters for the finite element method. The additional
  94 --   type parameter @n@ should be a type-level representation of the
  95 --   largest element in @max_degrees@. It needs to be known statically
  96 --   for the dimensions of the pointer matrix. The parameter @l@ is
  97 --   the number of global basis functions. It's equal to the number of
  98 --   /internal/ mesh nodes (i.e. m-1), plus the sum of (p_i - 1) for
  99 --   each p_i in max_degrees.
 100 data Params m n l a =
 101   Params {
 102     -- | A partition of the domain.
 103     mesh :: Col m (Interval a),
 104
 105     -- | A list of the highest-degree polynomial that we will use over
 106     --   each interval in our mesh.
 107     max_degrees :: Col m Int }
 108
 109
 110 -- | The pointer matrix mapping local basis elements to global
 111 --   ones. The (i,j)th entry of the matrix contains the index of the
 112 --   global basis function corresponding to the jth local basis
 113 --   function over element i.
 114 --
 115 --   TODO: This works for Dirichlet boundary conditions only.
 116 --
 117 --   Note that the number of columns in the matrix depends on the
 118 --
 119 --   Examples:
 120 --
 121 --   >>> import Data.Vector.Fixed ( N5, N6 )
 122 --   >>> import Linear.Matrix ( Col5, fromList )
 123 --   >>> import Naturals ( N19 )
 124 --
 125 --   >>> let p = fromList [[3],[3],[5],[4],[5]] :: Col5 Int
 126 --   >>> let mesh = undefined :: Col5 (Int,Int)
 127 --   >>> let params = Params mesh p :: Params N5 N5 N19 Int
 128 --   >>> let row1 = [0,1,5,6,0,0] :: [Int]
 129 --   >>> let row2 = [1,2,7,8,0,0] :: [Int]
 130 --   >>> let row3 = [2,3,9,10,11,12] :: [Int]
 131 --   >>> let row4 = [3,4,13,14,15,0] :: [Int]
 132 --   >>> let row5 = [4,0,16,17,18,19] :: [Int]
 133 --   >>> let expected = fromList [row1,row2,row3,row4,row5] :: Mat N5 N6 Int
 134 --   >>> pointer params == expected
 135 --   True
 136 --
 137 pointer :: (Arity m, Arity n, Arity l) => Params m n l a -> Mat m (S n) Int
 138 pointer params =
 139   construct lambda
 140   where
 141     -- | The number of polynomials in the local basis for element i.
 142     numpolys :: Int -> Int
 143     numpolys i = ((max_degrees params) !!! (i,0)) + 1
 144
 145     lambda i j
 146
 147       -- The first two (linear) basis functions are easy to do.
 148       | j < 2 = if i + j >= (nrows $ mesh params) then 0 else i + j
 149
 150       -- No local basis function exists for this j.
 151       | j >= (numpolys i) = 0
 152
 153       -- The first row we handle as a special case.
 154       | i == 0 = (nrows $ mesh params) + j - 2
 155
 156       -- The first entry for this row not corresponding to one of the
 157       -- linear polynomials (which we did in the first guard). This
 158       -- grabs the biggest number in the previous row and begins counting
 159       -- from there.
 160       | j == 2 = (lambda (i-1) (numpolys (i-1) - 1)) + 1
 161
 162       -- If none of the other cases apply, we can compute our value by
 163       -- looking to the left in the matrix.
 164       | otherwise = (lambda i (j-1)) + 1
 165
 166
 167
 168 -- | An affine transform taking the interval [x1,x2] into [-1,1].
 169 --
 170 --   Examples:
 171 --
 172 --   >>> let phi = affine (-6,7)
 173 --   >>> phi (-6)
 174 --   -1.0
 175 --   >>>  phi 7
 176 --   1.0
 177 --
 178 affine :: Field.C a => (a,a) -> (a -> a)
 179 affine (x1,x2) x = (fromInteger 2)*(x - x1)/(x2 - x1) - (fromInteger 1)
 180
 181 -- | The inverse of 'affine'. It should send [-1,1] into [x1,x2].
 182 --
 183 --   Examples:
 184 --
 185 --   >>> let phi = affine_inv (-6,7)
 186 --   >>> phi (-1)
 187 --   -6.0
 188 --   >>> phi 1
 189 --   7.0
 190 --
 191 affine_inv :: Field.C a => (a,a) -> (a -> a)
 192 affine_inv (x1,x2) x =
 193   x*(x2 - x1)/two + (x1 + x2)/two
 194   where
 195     two = fromInteger 2
 196
 197
 198 -- * Load vector
 199
 200 -- | Normalized integrals of orthogonal basis functions over
 201 --   [-1,1]. The test case below comes from Sage where the
 202 --   orthogonality of the polynomials' derivatives can easily be
 203 --   tested.
 204 --
 205 --   Examples:
 206 --
 207 --   >>> import qualified Algebra.Absolute as Absolute ( abs )
 208 --
 209 --   >>> let expected = 2.99624907925257
 210 --   >>> let actual = big_N 4 1.5 :: Double
 211 --   >>> Absolute.abs (actual - expected) < 1e-12
 212 --   True
 213 --
 214 big_N :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a
 215 big_N k x
 216   | k < 0 = error "requested a negative basis function"
 217   | k == 0 = (one - x) / (fromInteger 2)
 218   | k == 1 = (one + x) / (fromInteger 2)
 219   | otherwise =
 220       coeff * ( legendre k x - legendre (k-2) x )
 221       where
 222         two = fromInteger 2
 223         four = fromInteger 4
 224         coeff = one / (sqrt (four*(fromInteger k) - two)) :: a
 225
 226
 227 -- | A matrix containing 'big_N' functions indexed by their
 228 --   element/number. Each row in the matrix represents a finite element
 229 --   (i.e. an interval in the mesh). Within row @i@, column @j@ contains
 230 --   the @j@th 'big_N' basis function.
 231 --
 232 --   Any given 'big_N' will probably wind up in this matrix multiple
 233 --   times; the basis functions don't change depending on the
 234 --   interval. Only the /number/ of basis functions does. Computing
 235 --   them this way allows us to easily construct a lookup \"table\" of
 236 --   the proper dimensions.
 237 --
 238 --   The second example below relies on the fact that @big_N 3@ and
 239 --   @big_N 6@ expand to Legendre polynomials (2,4) and (5,7)
 240 --   respectively and so should be orthogonal over [-1,1].
 241 --
 242 --   Examples:
 243 --
 244 --   >>> import Data.Vector.Fixed ( N5, N6 )
 245 --   >>> import Integration.Gaussian ( gaussian )
 246 --   >>> import Linear.Matrix ( Col5, fromList )
 247 --   >>> import Naturals ( N19 )
 248 --
 249 --   >>> let p = fromList [[3],[3],[5],[4],[5]] :: Col5 Int
 250 --   >>> let mesh = undefined :: Col5 (Double,Double)
 251 --   >>> let params = Params mesh p :: Params N5 N5 N19 Double
 252 --   >>> let big_ns = big_N_matrix :: Mat N5 N6 (Double -> Double)
 253 --   >>> let n1 = big_ns !!! (1,0)
 254 --   >>> let n4 = big_ns !!! (4,0)
 255 --   >>> n1 1.5 == n4 1.5
 256 --   True
 257 --   >>> let n1 = big_ns !!! (1,3)
 258 --   >>> let n2 = big_ns !!! (2,4)
 259 --   >>> gaussian (\x -> (n1 x) * (n2 x)) < 1e-12
 260 --   True
 261 --
 262 big_N_matrix :: (Arity m, Arity n, Algebraic.C a, RealField.C a)
 263              => Mat m n (a -> a)
 264 big_N_matrix =
 265   construct lambda
 266   where
 267     lambda _ j x = big_N (toInteger j) x
 268
 269
 270 -- | The matrix of (N_i * N_j) functions used in the integrand of
 271 --   the mass matrices.
 272 big_Ns_matrix :: (Arity m, Arity n, Algebraic.C a, RealField.C a)
 273                => Mat m n (a -> a)
 274 big_Ns_matrix =
 275   construct lambda
 276   where
 277     lambda i j x = (big_N (toInteger i) x) * (big_N (toInteger j) x)
 278
 279
 280 -- | Compute the global load vector F.
 281 --
 282 --   Examples:
 283 --
 284 --   >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
 285 --   >>> import Naturals ( N3, N4, N7 )
 286 --
 287 --   >>> let big_A = const (1::Double)
 288 --   >>> let c x = sin x
 289 --   >>> let f x = x*(sin x)
 290 --   >>> let bdy = Left (Dirichlet (0,1::Double))
 291 --   >>> let pde = PDE big_A c f bdy
 292 --
 293 --   >>> let i1 = (0.0,1/3)
 294 --   >>> let i2 = (1/3,2/3)
 295 --   >>> let i3 = (2/3,4/5)
 296 --   >>> let i4 = (4/5,1.0)
 297 --   >>> let mesh = fromList [[i1], [i2], [i3], [i4]] :: Col4 (Double,Double)
 298 --   >>> let pvec = fromList [[2],[3],[2],[1]] :: Col4 Int
 299 --   >>> let params = Params mesh pvec :: Params N4 N3 N7 Double
 300 --
 301 --   >>> let f1 = [0.0418]
 302 --   >>> let f2 = [0.0805]
 303 --   >>> let f3 = [0.1007]
 304 --   >>> let f4 = [-0.0045]
 305 --   >>> let f5 = [-0.0332]
 306 --   >>> let f6 = [-0.0054]
 307 --   >>> let f7 = [-0.0267]
 308 --   >>> let expected = fromList [f1,f2,f3,f4,f5,f6,f7] :: Col N7 Double
 309 --   >>> let actual = big_F pde params
 310 --   >>> frobenius_norm (actual - expected) < 1e-4
 311 --   True
 312 --
 313 big_F :: forall m n l a.
 314          (Arity l, Arity m, Arity n,
 315           Algebraic.C a, RealField.C a, ToRational.C a)
 316       => PDE a
 317       -> Params m n l a
 318       -> Col l a
 319 big_F pde params =
 320   ifoldl2 accum zero (big_N_matrix :: Mat m (S n) (a -> a))
 321   where
 322     accum :: Int -> Int -> Col l a -> (a -> a) -> Col l a
 323     accum i j prev_F this_N =
 324       prev_F + this_F
 325       where
 326         two = fromInteger 2
 327         (x1,x2) = (mesh params) !!! (i,0)
 328         q = affine_inv (x1,x2)
 329         integrand x = ((f pde) (q x)) * (this_N x)
 330
 331         -- The pointer matrix numbers from 1 so subtract one here to
 332         -- get the right index.
 333         k = ((pointer params) !!! (i,j)) - 1
 334         integral = (gaussian integrand)*(x2 - x1) / two
 335         this_F = set_idx zero (k,0) integral
 336
 337
 338 -- * Stiffness matrix
 339
 340 -- | Derivatives of the 'big_N's, that is, orthogonal basis functions
 341 --   over [-1,1]. The test case below comes from Sage where the
 342 --   orthogonality of the polynomials' derivatives can easily be
 343 --   tested. The indices are shifted by one so that k=0 is the first
 344 --   basis function.
 345 --
 346 --   Examples:
 347 --
 348 --   >>> import qualified Algebra.Absolute as Absolute ( abs )
 349 --
 350 --   >>> let expected = 11.5757525403319
 351 --   >>> let actual = big_N' 3 1.5 :: Double
 352 --   >>> Absolute.abs (actual - expected) < 1e-10
 353 --   True
 354 --
 355 big_N' :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a
 356 big_N' k x
 357   | k < 0 = error "requested a negative basis function"
 358   | k == 0 = negate ( one / (fromInteger 2))
 359   | k == 1 = one / (fromInteger 2)
 360   | otherwise = coeff * ( legendre k x )
 361       where
 362         two = fromInteger 2
 363         coeff = sqrt ((two*(fromInteger k) + one) / two) :: a
 364
 365
 366 -- | The matrix of (N_i' * N_j') functions used in the integrand of
 367 --   the stiffness matrix.
 368 big_N's_matrix :: (Arity m, Arity n, Algebraic.C a, RealField.C a)
 369                => Mat m n (a -> a)
 370 big_N's_matrix =
 371   construct lambda
 372   where
 373     lambda i j x = (big_N' (toInteger i) x) * (big_N' (toInteger j) x)
 374
 375
 376 big_K_elem :: forall m n l a b.
 377          (Arity l, Arity m, Arity n,
 378           Algebraic.C a, RealField.C a, ToRational.C a)
 379       => PDE a
 380       -> Params m n l a
 381       -> Int
 382       -> Int
 383       -> Mat l l a
 384       -> b
 385       -> Mat l l a
 386 big_K_elem pde params _ k cur_K _ =
 387   ifoldl2 accum cur_K (big_N's_matrix :: Mat m (S n) (a -> a))
 388   where
 389     accum :: Int -> Int -> Mat l l a -> (a -> a) -> Mat l l a
 390     accum i j prev_K these_N's =
 391       prev_K + this_K
 392       where
 393         two = fromInteger 2
 394         (x1,x2) = (mesh params) !!! (k,0)
 395         q = affine_inv (x1,x2)
 396         integrand x = ((big_A pde) (q x)) * (these_N's x)
 397         -- The pointer matrix numbers from 1 so subtract one here to
 398         -- get the right index.
 399         row_idx = ((pointer params) !!! (k,i)) - 1
 400         col_idx = ((pointer params) !!! (k,j)) - 1
 401         integral = (two/(x2 - x1))* (gaussian integrand)
 402         this_K = set_idx zero (row_idx, col_idx) integral
 403
 404
 405
 406 -- | Compute the \"big K\" stiffness matrix. There are three
 407 --   parameters needed for K, namely i,j,k so a fold over a matrix will
 408 --   not do. This little gimmick simulates a three-index fold by doing a
 409 --   two-index fold over a row of the proper dimensions.
 410 --
 411 --   Examples:
 412 --
 413 --   >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
 414 --   >>> import Naturals ( N3, N4, N7 )
 415 --
 416 --   >>> let big_A = const (1::Double)
 417 --   >>> let c x = sin x
 418 --   >>> let f x = x*(sin x)
 419 --   >>> let bdy = Left (Dirichlet (0,1::Double))
 420 --   >>> let pde = PDE big_A c f bdy
 421 --
 422 --   >>> let i1 = (0.0,1/3)
 423 --   >>> let i2 = (1/3,2/3)
 424 --   >>> let i3 = (2/3,4/5)
 425 --   >>> let i4 = (4/5,1.0)
 426 --   >>> let mesh = fromList [[i1], [i2], [i3], [i4]] :: Col4 (Double,Double)
 427 --   >>> let pvec = fromList [[2],[3],[2],[1]] :: Col4 Int
 428 --   >>> let params = Params mesh pvec :: Params N4 N3 N7 Double
 429 --
 430 --   >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
 431 --   >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
 432 --   >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
 433 --   >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
 434 --   >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
 435 --   >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
 436 --   >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
 437 --   >>> let expected = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
 438 --   >>> let actual = big_K pde params
 439 --   >>> frobenius_norm (actual - expected) < 1e-10
 440 --   True
 441 --
 442 big_K :: forall m n l a.
 443          (Arity l, Arity m, Arity n,
 444           Algebraic.C a, RealField.C a, ToRational.C a)
 445       => PDE a
 446       -> Params m n l a
 447       -> Mat l l a
 448 big_K pde params =
 449   ifoldl2 (big_K_elem pde params) zero col_idxs
 450   where
 451     col_idxs = fromList [map fromInteger [0..]] :: Row m a
 452
 453
 454 -- * Mass matrix
 455
 456 big_M_elem :: forall m n l a b.
 457          (Arity l, Arity m, Arity n,
 458           Algebraic.C a, RealField.C a, ToRational.C a)
 459       => PDE a
 460       -> Params m n l a
 461       -> Int
 462       -> Int
 463       -> Mat l l a
 464       -> b
 465       -> Mat l l a
 466 big_M_elem pde params _ k cur_M _ =
 467   ifoldl2 accum cur_M (big_Ns_matrix :: Mat m (S n) (a -> a))
 468   where
 469     accum :: Int -> Int -> Mat l l a -> (a -> a) -> Mat l l a
 470     accum i j prev_M these_Ns =
 471       prev_M + this_M
 472       where
 473         two = fromInteger 2
 474         (x1,x2) = (mesh params) !!! (k,0)
 475         q = affine_inv (x1,x2)
 476         integrand x = ((c pde) (q x)) * (these_Ns x)
 477         -- The pointer matrix numbers from 1 so subtract one here to
 478         -- get the right index.
 479         row_idx = ((pointer params) !!! (k,i)) - 1
 480         col_idx = ((pointer params) !!! (k,j)) - 1
 481         integral = (x2 - x1)*(gaussian integrand) / two
 482         this_M = set_idx zero (row_idx, col_idx) integral
 483
 484
 485 -- | Compute the \"big M\" mass matrix. There are three
 486 --   parameters needed for M, namely i,j,k so a fold over a matrix will
 487 --   not do. This little gimmick simulates a three-index fold by doing a
 488 --   two-index fold over a row of the proper dimensions.
 489 --
 490 --   Examples:
 491 --
 492 --   >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
 493 --   >>> import Naturals ( N3, N4, N7 )
 494 --
 495 --   >>> let big_A = const (1::Double)
 496 --   >>> let c x = sin x
 497 --   >>> let f x = x*(sin x)
 498 --   >>> let bdy = Left (Dirichlet (0,1::Double))
 499 --   >>> let pde = PDE big_A c f bdy
 500 --
 501 --   >>> let i1 = (0.0,1/3)
 502 --   >>> let i2 = (1/3,2/3)
 503 --   >>> let i3 = (2/3,4/5)
 504 --   >>> let i4 = (4/5,1.0)
 505 --   >>> let mesh = fromList [[i1], [i2], [i3], [i4]] :: Col4 (Double,Double)
 506 --   >>> let pvec = fromList [[2],[3],[2],[1]] :: Col4 Int
 507 --   >>> let params = Params mesh pvec :: Params N4 N3 N7 Double
 508 --
 509 --   >>> let m1 = [0.0723,0.0266,0,-0.0135,-0.0305,0.0058,0] :: [Double]
 510 --   >>> let m2 = [0.0266,0.0897,0.0149,0,-0.0345,-0.0109,-0.0179] :: [Double]
 511 --   >>> let m3 = [0,0.0149,0.0809,0,0,0,-0.0185] :: [Double]
 512 --   >>> let m4 = [-0.0135,0,0,0.0110,0,0,0] :: [Double]
 513 --   >>> let m5 = [-0.0305,-0.0345,0,0,0.0319,0.0018,0] :: [Double]
 514 --   >>> let m6 = [0.0058,-0.0109,0,0,0.0018,0.0076,0] :: [Double]
 515 --   >>> let m7 = [0,-0.0179,-0.0185,0,0,0,0.0178] :: [Double]
 516 --
 517 --   >>> let expected = fromList [m1,m2,m3,m4,m5,m6,m7] :: Mat N7 N7 Double
 518 --   >>> let actual = big_M pde params
 519 --   >>> frobenius_norm (actual - expected) < 1e-3
 520 --   True
 521 --
 522 big_M :: forall m n l a.
 523          (Arity l, Arity m, Arity n,
 524           Algebraic.C a, RealField.C a, ToRational.C a)
 525       => PDE a
 526       -> Params m n l a
 527       -> Mat l l a
 528 big_M pde params =
 529   ifoldl2 (big_M_elem pde params) zero col_idxs
 530   where
 531     col_idxs = fromList [map fromInteger [0..]] :: Row m a
 532
 533
 534
 535 -- | Determine the coefficient vector @x@ from the system @(K + M)x = F@.
 536 --
 537 --   Examples:
 538 --
 539 --   >>> import Linear.Matrix ( Col4, Col7, frobenius_norm, fromList )
 540 --   >>> import Naturals ( N3, N4, N7 )
 541 --
 542 --   >>> let big_A = const (1::Double)
 543 --   >>> let c x = sin x
 544 --   >>> let f x = x*(sin x)
 545 --   >>> let bdy = Left (Dirichlet (0,1::Double))
 546 --   >>> let pde = PDE big_A c f bdy
 547 --
 548 --   >>> let i1 = (0.0,1/3)
 549 --   >>> let i2 = (1/3,2/3)
 550 --   >>> let i3 = (2/3,4/5)
 551 --   >>> let i4 = (4/5,1.0)
 552 --   >>> let mesh = fromList [[i1], [i2], [i3], [i4]] :: Col4 (Double,Double)
 553 --   >>> let pvec = fromList [[2],[3],[2],[1]] :: Col4 Int
 554 --   >>> let params = Params mesh pvec :: Params N4 N3 N7 Double
 555 --
 556 --   >>> let c1 = [0.02366220347687] :: [Double]
 557 --   >>> let c2 = [0.03431630082636] :: [Double]
 558 --   >>> let c3 = [0.02841800893264] :: [Double]
 559 --   >>> let c4 = [-0.00069489654996] :: [Double]
 560 --   >>> let c5 = [-0.00518637005151] :: [Double]
 561 --   >>> let c6 = [-0.00085028505337] :: [Double]
 562 --   >>> let c7 = [-0.00170478210110] :: [Double]
 563 --   >>> let expected = fromList [c1,c2,c3,c4,c5,c6,c7] :: Col7 Double
 564 --   >>> let actual = coefficients pde params
 565 --   >>> frobenius_norm (actual - expected) < 1e-8
 566 --   True
 567 --
 568 coefficients :: forall m n l a.
 569                 (Arity m, Arity n, Arity l,
 570                  Algebraic.C a, Eq a, RealField.C a, ToRational.C a)
 571              => PDE a
 572              -> Params m n (S l) a
 573              -> Col (S l) a
 574 coefficients pde params =
 575   solve_positive_definite matrix b
 576   where
 577     matrix = (big_K pde params) + (big_M pde params)
 578     b      = big_F pde params
 579
 580
 581 solution :: forall m n l a.
 582             (Arity m, Arity n, Arity l,
 583             Algebraic.C a, Eq a, RealField.C a, ToRational.C a, Show a)
 584          => PDE a
 585          -> Params m n (S l) a
 586          -> Piecewise a
 587 solution pde params =
 588   from_intervals $ map head $ toList $ solved_column
 589   where
 590     global_coeffs :: Col (S l) a
 591     global_coeffs = coefficients pde params
 592
 593     ptr :: Mat m (S n) Int
 594     ptr = pointer params
 595
 596     -- Each mesh element has an associated row in the pointer
 597     -- matrix. Stick them together.
 598     mesh_with_ptr_rows :: Col m (Interval a, Row (S n) Int)
 599     mesh_with_ptr_rows = zip2 (mesh params) (rows2 ptr)
 600
 601     make_local_coeffs :: (Interval a, Row (S n) Int) -> Row (S n) a
 602     make_local_coeffs (interval, ptr_row) =
 603       construct lambda
 604         where
 605           lambda _ j = if (ptr_row !!! (0,j)) == zero
 606                        then zero
 607                        else global_coeffs !!! ((ptr_row !!! (0,j)) - 1, 0)
 608
 609     -- Create a column vector for each mesh element containing the global
 610     -- coefficients corresponding to that element.
 611     local_coeffs :: Col m (Row (S n) a)
 612     local_coeffs = map2 make_local_coeffs mesh_with_ptr_rows
 613
 614     global_basis_functions :: Col (S n) (a -> a)
 615     global_basis_functions =
 616       construct lambda
 617         where lambda i _ = big_N (toInteger i)
 618
 619     mesh_with_coeffs :: Col m (Interval a, Row (S n) a)
 620     mesh_with_coeffs = zip2 (mesh params) local_coeffs
 621
 622     solved_column :: Col m (Interval a, (a -> a))
 623     solved_column = map2 solve_piece $ mesh_with_coeffs
 624
 625     solve_piece :: (Interval a, Row (S n) a) -> (Interval a, (a -> a))
 626     solve_piece (interval, coeffs_row) = (interval, f)
 627       where
 628         coeffs_col = transpose coeffs_row
 629
 630         f x = element_sum2 $ zipwith2 combine coeffs_col global_basis_functions
 631           where
 632             xi = (affine interval) x
 633             combine ci ni = ci*(ni xi)
 634