Start page of documentation

Download the code

Documentation: PDF, PS

Additional equations for jumps

At first, jumps are expressed in terms of known quantities, like boundary conditions, and one sided derivatives. Then, for jumps at arbitrary intersection point we can write equations

• Dirichlet b.c.

  $$\begin{array}{rcccl} && [u] &=& -u_D \partial_x u^- & + & ... ...0 \partial_{xy} u^- & + &[\partial_{xy} u] &=& 0 \end{array}$$ (8)

  
  

• Neumann b.c.

  $$\begin{array}{rcccl} u^- & + &[u] &=& 0 t_1^2 \partial_x u... ...0 \partial_{xy} u^- & +& [\partial_{xy} u] &=& 0 \end{array}$$ (9)

  
  

Note that the one sided derivatives depend on the unknown solution. Our goal is to approximate them using a least squares fit of a quadratic polynomial.

The first step is to select the stencil for least squares interpolation. We follow the methodology from [3]. It is done by selecting some radius $k\in\mathbb{N}$ and then taking all grid points of $\Omega^-$ inside circle with radius $k \sqrt{h_x h_y}$. (Line 73 in DIFFOP/d_matrix.m.) With $l$ we denote the stencil cardinality.

For each stencil point we try to fit a quadratic polynomial:

$$u(x_i,y_j)\approx p(x_i,y_j) :=p_1 + p_2 h_i + p_3 k_j + p_4 h_i^2 + p_5 k_j^2 + p_6 h_i k_j ,$$

where $h_i:=x_i-x_\alpha$ and $k_j:=y_j-y_\alpha$. With this a weighted least squares problem is obtained:

$$\sum_{(x_i,y_j)\in\mbox{ stencil}} w_{i,j}^2 \left( p(x_i,y_j) - u_{i,j} \right)^2 \xrightarrow[P]{} \min$$

where $P:=(p_1, p_2, p_3, p_4, p_5, p_6)^\top$. In the code the following weights are used:

$$w_{i,j} = \left(1 + d(x_i,y_j)/h_x \right)^{-1} , d(x_i,y_j):=\sqrt{h_i^2 + k_j^2} .$$

For a given grid function $U$, coefficients of the polynomial are given by

$$P = (M^\top W^2 M)^{-1}(M^\top W^2 R U) ,$$

where

$$W:= \mathrm{diag}(w_{i_1,j_1}, w_{i_2,j_2}, \ldots, w_{i_l,j_l})$$

and

$$M:= \begin{pmatrix} \vdots & \vdots & \vdots & \vdots & \vdots & ... ... h_i k_j \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \end{pmatrix}$$

(in code this is the matrix MM in DIFFOP/d_matrix.m.)

Matrix $R$ is the restriction operator, $RU$ restricts the vector $U$ defined at all grid points to its values only at the stencil points. $R$ is possible to write explicitely, however, we have found that this costs enormous computational time and because of this reason in code the matrix $W^2 R$ is directly computed (matrix W2R in DIFFOP/d_matrix.m).

Note that

$$u(x_\alpha,y_\alpha) \approx p_1$$

$$\partial_x u(x_\alpha,y_\alpha) \approx p_2$$

$$\partial_y u(x_\alpha,y_\alpha) \approx p_3$$

$$\partial_{xx} u(x_\alpha,y_\alpha) \approx 2p_4$$

$$\partial_{yy} u(x_\alpha,y_\alpha) \approx 2p_5$$

$$\partial_{xy} u(x_\alpha,y_\alpha) \approx p_6$$

and this allows us to write the corresponding one-sided derivatives as

$$(u^-, \partial_x u^-, \partial_y u^-, \partial_{xx} u^-, \partial_{yy} u^-, \partial_{xy} u^-)^\top = BU ,$$

where

$$B:=S (M\top W^2 M)^{-1} (M^\top W^2 R) .$$

In DIFFOP/d_matrix.m to matrix $B$ corresponds the variable B.

Using (8) we obtain in the case of Dirichlet boundary condition

$$\begin{array}{cccccc} 0\cdot B(1,:) &U &+& [u] &=& -u_D \\ ... ...mbox{These form vector } \tilde F}{\underbrace{ 0}} \end{array}$$

Using (9) we get for Neumann boundary

$$\begin{array}{cccccc} 1\cdot B(1,:)& U &+& [u] &=& 0 \\ \le... ...mbox{These form vector } \tilde F}{\underbrace{ 0}} \end{array}$$

In code this corresponds to lines 95--115 in DIFFOP/d_matrix.m.

Thus, the system for jumps can be written as

$$\mathbf{D}U + J = \tilde F .$$ (10)

In code, the matrix $\mathbf{D}$ (stored in variable D) and the right hand side vector $\tilde F$ (stored in variable Ft) consist each of two parts, corresponding to Dirichlet and Neumann intersection points:

$$\mathbf{D} = \begin{pmatrix}\mathbf{D}_D \mathbf{D}_D \end{pma... ... , \ \tilde F = \begin{pmatrix}\tilde F_D \tilde F_N \end{pmatrix} .$$

(ASSEMBLE/run_ejiim.m.)

Note: matrix $\mathbf{D}$ and vector $\tilde F$ depend on boundary conditions!