Poisson Equation

\[\begin{split}\begin{cases} -\nabla^2 u = f(x) &\text{ in } \Omega \\ u = 0 & \text{ on } \partial\Omega \end{cases}\end{split}\]
mozart.poisson.solve.DVandermondeM1D(degree, r)[source]

Initialize the derivative of modal basis (i) at (r) at order degree

Paramters
  • degree (int32) : Polynomial degree
  • r (float64 array) : points
Returns
  • DVr (float64 array) : Differentiate Vandermonde matrix
Example
>>> N = 3
>>> from mozart.poisson.solve import VandermondeM1D
>>> r = np.linspace(-1,1,N+1)
>>> DVr = DVandermondeM1D(N,r)
>>> print(DVr)
array([[  0.        ,   1.22474487,  -4.74341649,  11.22497216],
[  0.        ,   1.22474487,  -1.58113883,  -1.24721913],
[  0.        ,   1.22474487,   1.58113883,  -1.24721913],
[  0.        ,   1.22474487,   4.74341649,  11.22497216]])
mozart.poisson.solve.Dmatrix1D(degree, r, V)[source]

Initialize the derivative of modal basis (i) at (r) at order degree

Paramters
  • degree (int32) : Polynomial degree
  • r (float64 array) : points
Returns
  • DVr (float64 array) : Differentiate Vandermonde matrix
Example
>>> N = 3
>>> from mozart.poisson.solve import Dmatrix1D
>>> r = np.linspace(-1,1,N+1)
>>> Dr = Dmatirx1D(N,r)
>>> print(Dr)
array([[-2.75,  4.5 , -2.25,  0.5 ],
[-0.5 , -0.75,  1.5 , -0.25],
[ 0.25, -1.5 ,  0.75,  0.5 ],
[-0.5 ,  2.25, -4.5 ,  2.75]])
mozart.poisson.solve.DnJacobiP(x, alpha=0, beta=0, degree=0)[source]

Evaluate the derivative of the normalized Jacobi polynomial of type alpha, beta > -1 at point x for order n

Paramters
  • x (float64 array) : variable x
  • alpha (int32) : superscript alpha of normalized Jacobi polynomial
  • beta (int32) : superscript beta of normalized Jacobi polynomial
  • degree (int32) : Polynomial degree
Returns
  • dP (float64 array) : the value of the derivative of the normalized Jacobi polynomial at x according to alpha, beta, degree
Example
>>> N = 2
>>> x = np.array([-1, 0, 1])
>>> from mozart.poisson.solve import DnJacobiP
>>> dP = DnJacobiP(x,0,0,N)
>>> print(dP)
array([-4.74341649,  0.        ,  4.74341649])
mozart.poisson.solve.VandermondeM1D(degree, r)[source]

Initialize the 1D Vandermonde matrix, \(V_{i,j} = \phi_j(r_i)\)

Paramters
  • degree (int32) : Polynomial degree
  • r (float64 array) : points
Returns
  • V1D (float64 array) : 1D Vandermonde matrix
Example
>>> N = 3
>>> from mozart.poisson.solve import VandermondeM1D
>>> r = np.linspace(-1,1,N+1)
>>> V1D = VandermondeM1D(N,r)
>>> print(V1D)
array([[ 0.70710678, -1.22474487,  1.58113883, -1.87082869],
[ 0.70710678, -0.40824829, -0.52704628,  0.76218947],
[ 0.70710678,  0.40824829, -0.52704628, -0.76218947],
[ 0.70710678,  1.22474487,  1.58113883,  1.87082869]])
mozart.poisson.solve.computeError_one_dim(c4n, n4e, ind4e, exact_u, exact_ux, approx_u, degree, degree_i)[source]

Computes L^2-error and semi H^1-error between exact solution and approximate solution.

Parameters
  • c4n (float64 array) : coordinates for nodes
  • n4e (int32 array) : nodes for elements
  • ind4e (int32 array) : indices for elements
  • exact_u (lambda) : exact solution
  • exact_ux (lambda) : derivative of exact solution
  • approx_u (float64 array) : approximate solution
  • degree (int32) : Polynomial degree
  • degree_i (int32) : Polynomial degree for interpolation
Returns
  • L2error (float64) : L^2 error between exact solution and approximate solution.
  • sH1error (float64) : semi H^1 error between exact solution and approximate solution.
Example
>>> N = 2
>>> from mozart.mesh.rectangle import interval 
>>> c4n, n4e, n4db, ind4e = interval(0, 1, 4, 2)
>>> f = lambda x: np.pi ** 2 * np.sin(np.pi * x)
>>> u_D = lambda x: np.zeros_like(x)
>>> from mozart.poisson.solve import one_dim_p
>>> x = one_dim_p(c4n, n4e, n4db, ind4e, f, u_D, N)
>>> from mozart.poisson.solution import computeError_one_dim
>>> exact_u = lambda x: np.sin(np.pi * x)
>>> exact_ux = lambda x: np.pi * np.cos(np.pi * x)
>>> L2error, sH1error = computeError_one_dim(c4n, n4e, ind4e, exact_u, exact_ux, x, N, N+3)
>>> L2error
0.0020225729623142077
>>> sH1error
0.05062779815975444
mozart.poisson.solve.getMatrix1D(degree)[source]

Get FEM matrices on the reference domain I = [-1, 1]

Paramters
  • degree (int32) : degree of polynomial
Returns
  • M_R (float64 array) : Mass matrix on the reference domain
  • S_R (float64 array) : Stiffness matrix on the reference domain
  • D_R (float64 array) : Differentiation matrix on the reference domain
mozart.poisson.solve.nJacobiGL(alpha=0, beta=0, degree=0)[source]

Compute the degree-th order Gauss Lobatto quadrature points x associated with the nomalized Jacobi polynomial of type \(\alpha, \beta > -1\)

Paramters
  • alpha (int32) : superscript alpha of normalized Jacobi polynomial
  • beta (int32) : superscript beta of normalized Jacobi polynomial
  • degree (int32) : Polynomial degree
Returns
  • x (float64 array) : Gauss Lobatto quadrature points
Example
>>> N = 3
>>> from mozart.poisson.solve import nJacobiGL
>>> x = nJacobiGQ(0,0,N)
>>> print(x)
array([-1.       , -0.4472136,  0.4472136,  1.       ])
mozart.poisson.solve.nJacobiGQ(alpha=0, beta=0, degree=0)[source]

Compute the degree-th order Gauss quadrature points x and weights w associated with the nomalized Jacobi polynomial of type alpha, beta > -1

Paramters
  • alpha (int32) : superscript alpha of normalized Jacobi polynomial
  • beta (int32) : superscript beta of normalized Jacobi polynomial
  • degree (int32) : Polynomial degree
Returns
  • x (float64 array) : Gauss quadrature points
  • w (float64 array) : Gauss quadrature weights
Example
>>> N = 2
>>> from mozart.poisson.solve import nJacobiGQ
>>> x, w = nJacobiGQ(0,0,N)
>>> print(x)
array([ -7.74596669e-01,  -4.78946310e-17,   7.74596669e-01])
>>> print(w)
array([ 0.55555556,  0.88888889,  0.55555556])
mozart.poisson.solve.nJacobiP(x, alpha=0, beta=0, degree=0)[source]

Evaluate normalized Jacobi polynomial of type alpha, beta > -1 at point x for order n (the special case of alpha = beta = 0, knows as the normalized Legendre polynomial)

Paramters
  • x (float64 array) : variable x
  • alpha (int32) : superscript alpha of normalized Jacobi polynomial
  • beta (int32) : superscript beta of normalized Jacobi polynomial
  • degree (int32) : Polynomial degree
Returns
  • P (float64 array) : the value of degree-th order normalized Jacobi polynomial at x
Example
>>> N = 2
>>> x = np.array([-1, 0, 1])
>>> from mozart.poisson.solve import nJacobiP
>>> P = nJacobiP(x,0,0,N)
>>> print(P)
array([ 1.58113883, -0.79056942,  1.58113883])
mozart.poisson.solve.one_dim(c4n, n4e, n4Db, f, u_D, degree=1)[source]

Computes the coordinates of nodes and elements.

Parameters
  • c4n (float64 array) : coordinates
  • n4e (int32 array) : nodes for elements
  • n4Db (int32 array) : Dirichlet boundary nodes
  • f (lambda) : source term
  • u_D (lambda) : Dirichlet boundary condition
  • degree (int32) : Polynomial degree
Returns
  • x (float64 array) : solution
Example
>>> N = 3
>>> c4n, n4e = unit_interval(N)
>>> n4Db = [0, N-1]
>>> f = lambda x: np.ones_like(x)
>>> u_D = lambda x: np.zeros_like(x)
>>> from mozart.poisson.solve import one_dim
>>> x = one_dim(c4n, n4e, n4Db, f, u_D)
>>> print(x)
array([ 0.   ,  0.125,  0.   ])
mozart.poisson.solve.one_dim_p(c4n, n4e, n4db, ind4e, f, u_D, degree)[source]

Computes the coordinates of nodes and elements.

Parameters
  • c4n (float64 array) : coordinates for nodes
  • n4e (int32 array) : nodes for elements
  • n4db (int32 array) : nodes for Dirichlet boundary
  • ind4e (int32 array) : indices for elements
  • f (lambda) : source term
  • u_D (lambda) : Dirichlet boundary condition
  • degree (int32) : Polynomial degree
Returns
  • x (float64 array) : solution
Example
>>> N = 2
>>> from mozart.mesh.rectangle import interval 
>>> c4n, n4e, n4db, ind4e = interval(0, 1, 4, 2)
>>> f = lambda x: np.ones_like(x)
>>> u_D = lambda x: np.zeros_like(x)
>>> from mozart.poisson.solve import one_dim_p
>>> x = one_dim_p(c4n, n4e, n4db, ind4e, f, u_D, N)
>>> x
array([ 0.       ,  0.0546875,  0.09375  ,  0.1171875,  0.125    ,
   0.1171875,  0.09375  ,  0.0546875,  0.       ])