2 {-# LANGUAGE ScopedTypeVariables #-}
4 -- | Finite Element Method (FEM) to solve the problem,
6 -- -(A(x)*u'(x))' + c(x)*u*(x) = f(x) on (a,b)
8 -- with one of the two types of boundary conditions,
10 -- * Dirichlet: u(a) = u(b) = 0
12 -- * A(a)u'(a) = alpha, A(b)u'(b) = beta
14 -- if c(x) = 0, then it is assumed that the boundary conditions are
17 -- The code creates a linear system,
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.
25 -- Warning: until PLU factorization is implemented, we can only
26 -- solve the resulting system if it's positive definite!
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 )
37 import qualified Prelude as P
39 import Integration.Gaussian ( gaussian )
40 import Linear.Matrix (
57 import Linear.System ( solve_positive_definite )
58 import Piecewise ( Piecewise(..), from_intervals )
59 import Polynomials.Orthogonal ( legendre )
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 }
65 -- | Neumann boundary conditions. @alpha@ specifies A(a)u'(b) and
66 -- @beta@ specifies A(b)u'(b).
68 Neumann { domain_neumann :: Interval a,
72 -- | Boundary conditions can be either Dirichlet or Neumann.
73 type BoundaryConditions a = Either (Dirichlet a) (Neumann a)
75 type Interval a = (a,a)
77 -- | The data for the PDE that we are attempting to solve.
87 -- | The boundary conditions. The domain also specifies the
89 bdy :: BoundaryConditions a }
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 =
102 -- | A partition of the domain.
103 mesh :: Col m (Interval a),
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 }
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.
115 -- TODO: This works for Dirichlet boundary conditions only.
117 -- Note that the number of columns in the matrix depends on the
121 -- >>> import Data.Vector.Fixed ( N5, N6 )
122 -- >>> import Linear.Matrix ( Col5, fromList )
123 -- >>> import Naturals ( N19 )
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
137 pointer :: (Arity m, Arity n, Arity l) => Params m n l a -> Mat m (S n) Int
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
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
150 -- No local basis function exists for this j.
151 | j >= (numpolys i) = 0
153 -- The first row we handle as a special case.
154 | i == 0 = (nrows $ mesh params) + j - 2
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
160 | j == 2 = (lambda (i-1) (numpolys (i-1) - 1)) + 1
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
168 -- | An affine transform taking the interval [x1,x2] into [-1,1].
172 -- >>> let phi = affine (-6,7)
178 affine :: Field.C a => (a,a) -> (a -> a)
179 affine (x1,x2) x = (fromInteger 2)*(x - x1)/(x2 - x1) - (fromInteger 1)
181 -- | The inverse of 'affine'. It should send [-1,1] into [x1,x2].
185 -- >>> let phi = affine_inv (-6,7)
191 affine_inv :: Field.C a => (a,a) -> (a -> a)
192 affine_inv (x1,x2) x =
193 x*(x2 - x1)/two + (x1 + x2)/two
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
207 -- >>> import qualified Algebra.Absolute as Absolute ( abs )
209 -- >>> let expected = 2.99624907925257
210 -- >>> let actual = big_N 4 1.5 :: Double
211 -- >>> Absolute.abs (actual - expected) < 1e-12
214 big_N :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a
216 | k < 0 = error "requested a negative basis function"
217 | k == 0 = (one - x) / (fromInteger 2)
218 | k == 1 = (one + x) / (fromInteger 2)
220 coeff * ( legendre k x - legendre (k-2) x )
224 coeff = one / (sqrt (four*(fromInteger k) - two)) :: a
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.
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.
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].
244 -- >>> import Data.Vector.Fixed ( N5, N6 )
245 -- >>> import Integration.Gaussian ( gaussian )
246 -- >>> import Linear.Matrix ( Col5, fromList )
247 -- >>> import Naturals ( N19 )
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
257 -- >>> let n1 = big_ns !!! (1,3)
258 -- >>> let n2 = big_ns !!! (2,4)
259 -- >>> gaussian (\x -> (n1 x) * (n2 x)) < 1e-12
262 big_N_matrix :: (Arity m, Arity n, Algebraic.C a, RealField.C a)
267 lambda _ j x = big_N (toInteger j) x
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)
277 lambda i j x = (big_N (toInteger i) x) * (big_N (toInteger j) x)
280 -- | Compute the global load vector F.
284 -- >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
285 -- >>> import Naturals ( N3, N4, N7 )
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
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
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
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)
320 ifoldl2 accum zero (big_N_matrix :: Mat m (S n) (a -> a))
322 accum :: Int -> Int -> Col l a -> (a -> a) -> Col l a
323 accum i j prev_F this_N =
327 (x1,x2) = (mesh params) !!! (i,0)
328 q = affine_inv (x1,x2)
329 integrand x = ((f pde) (q x)) * (this_N x)
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
338 -- * Stiffness matrix
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
348 -- >>> import qualified Algebra.Absolute as Absolute ( abs )
350 -- >>> let expected = 11.5757525403319
351 -- >>> let actual = big_N' 3 1.5 :: Double
352 -- >>> Absolute.abs (actual - expected) < 1e-10
355 big_N' :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a
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 )
363 coeff = sqrt ((two*(fromInteger k) + one) / two) :: a
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)
373 lambda i j x = (big_N' (toInteger i) x) * (big_N' (toInteger j) x)
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)
386 big_K_elem pde params _ k cur_K _ =
387 ifoldl2 accum cur_K (big_N's_matrix :: Mat m (S n) (a -> a))
389 accum :: Int -> Int -> Mat l l a -> (a -> a) -> Mat l l a
390 accum i j prev_K these_N's =
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
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.
413 -- >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
414 -- >>> import Naturals ( N3, N4, N7 )
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
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
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
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)
449 ifoldl2 (big_K_elem pde params) zero col_idxs
451 col_idxs = fromList [map fromInteger [0..]] :: Row m a
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)
466 big_M_elem pde params _ k cur_M _ =
467 ifoldl2 accum cur_M (big_Ns_matrix :: Mat m (S n) (a -> a))
469 accum :: Int -> Int -> Mat l l a -> (a -> a) -> Mat l l a
470 accum i j prev_M these_Ns =
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
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.
492 -- >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
493 -- >>> import Naturals ( N3, N4, N7 )
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
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
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]
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
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)
529 ifoldl2 (big_M_elem pde params) zero col_idxs
531 col_idxs = fromList [map fromInteger [0..]] :: Row m a
535 -- | Determine the coefficient vector @x@ from the system @(K + M)x = F@.
539 -- >>> import Linear.Matrix ( Col4, Col7, frobenius_norm, fromList )
540 -- >>> import Naturals ( N3, N4, N7 )
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
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
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
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)
572 -> Params m n (S l) a
574 coefficients pde params =
575 solve_positive_definite matrix b
577 matrix = (big_K pde params) + (big_M pde params)
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)
585 -> Params m n (S l) a
587 solution pde params =
588 from_intervals $ map head $ toList $ solved_column
590 global_coeffs :: Col (S l) a
591 global_coeffs = coefficients pde params
593 ptr :: Mat m (S n) Int
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)
601 make_local_coeffs :: (Interval a, Row (S n) Int) -> Row (S n) a
602 make_local_coeffs (interval, ptr_row) =
605 lambda _ j = if (ptr_row !!! (0,j)) == zero
607 else global_coeffs !!! ((ptr_row !!! (0,j)) - 1, 0)
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
614 global_basis_functions :: Col (S n) (a -> a)
615 global_basis_functions =
617 where lambda i _ = big_N (toInteger i)
619 mesh_with_coeffs :: Col m (Interval a, Row (S n) a)
620 mesh_with_coeffs = zip2 (mesh params) local_coeffs
622 solved_column :: Col m (Interval a, (a -> a))
623 solved_column = map2 solve_piece $ mesh_with_coeffs
625 solve_piece :: (Interval a, Row (S n) a) -> (Interval a, (a -> a))
626 solve_piece (interval, coeffs_row) = (interval, f)
628 coeffs_col = transpose coeffs_row
630 f x = element_sum2 $ zipwith2 combine coeffs_col global_basis_functions
632 xi = (affine interval) x
633 combine ci ni = ci*(ni xi)