# Points on the interface

We consider the problem

$\left\{\begin{matrix} k(x)\Delta u(x)=f(x) & \text{ in } \Omega\\ u=0 & \text{ in } \Gamma \end{matrix}\right.$

where $\Omega \subset \mathbb{R}^2$ open and bounded with smooth boundary $\Gamma=\partial{\Omega}$ and $k(x) \geq a, \forall x \in \Omega$.

We consider that $\Omega= \Omega_1 \cup \Omega_2$ and that $k$ is always a constant $k(x)=\left\{\begin{matrix} k_1 & \text{ in } \Omega_1 \\ k_2 & \text{ in } \Omega_2 \end{matrix}\right.$

$\Sigma$ is the boundary between $\Omega_1$ and $\Omega_2$, $\Sigma= \overline{\Omega_1} \cap \overline{\Omega_2}, \Gamma_j=\Gamma \cap \partial{\Omega_j}, j=1,2$.

We suppose that $\Omega=[-1,1]^2$ and $\Sigma=[-1,1] \times \{ y=0 \}, \Omega_1=\Omega \cap \{ y>0 \}$ and $\Omega_2=\Omega \cap \{ y<0 \}$.

We construct a grid $\Omega$ so that the interface $\Sigma$ coincides with the sides of the squares.

We write the problem as $AU=f$ in order to find the approximation of the solution.

For the approximations that correspond to the points that are over $\Sigma$ we pick $k1$ and for the approximations that correspond to the points that are under $\Sigma$ we pick $k2$.

What do we do for the approimations that correspond to the points that are on $\Sigma$ ?

EDIT:Finally, we don't want the approximation of the solution to be continuous. We may take $k=\frac{k_1+k_2}{2}$ for the points on the interface. So do we have to take cases for $y$ in order to give a value to $k$?

• what physical problem do you considering? Nov 15 '15 at 22:51
• Is your interface moving? Nov 17 '15 at 22:40
• No, it's the line $y=0$. @nluigi Nov 17 '15 at 23:15
• In that case perhaps my answer to this question might be useful; specifically check out the section Simulating a composite material. It is based on a Finite Volume approach in which the interface is located halfway between nodes. Nov 18 '15 at 14:04
• @nluigi Could you look at my question of the last comment if you have an idea? Nov 19 '15 at 15:20

Rewrite your problem as $$\begin{cases} \Delta u(x)=\frac{f(x)}{k_1}& \text{ in } \Omega_1\\ \Delta u(x)=\frac{f(x)}{k_2}& \text{ in } \Omega_2\\ u=0 & \text{ in } \Gamma \end{cases}$$ you have to add a condition on an interface, for example a jump condition, e.g. $$u^1(x)-u^2(x)=\alpha(x)\quad\text{and}\quad u^1_n(x)-u^2_n(y)=\beta(x) ,\quad x\in\Sigma$$ or smoothness (e.g. $\alpha=0=\beta$) of the solution and its normal derivative across the interface. In the above $u^i$ is the solution in the domain $\Omega_i$, $i\in\{1,2\}$.

One of the best ways to solve such problem would be Difference Potentials Method (based on good theory, need time to learn, worth the time you spend on it), another is Immersed Interface method(easy to understand, implementation is involved and changes with the problem, less accurate). You find a little more under the link provided.

Regularization approach

If the low order of accuracy is enough for you, you may want to try something very simple. You take a continuous approximation of step function and solve a little bit different, but continuous problem, e.g.: $$\begin{cases} \Delta u(x)=\frac{f(x)}{g_p(x)}& \text{ in } \Omega_1\cup \Omega_2\\ u=0 & \text{ in } \Gamma \end{cases}$$ where (you can choose another approximation of step function) $$g_p=\frac{k_2-k_1}{1+e^{-M(x-p)}} + k_1.$$ $M$ supposed to be a big number, the larger the $M$ the sharper the "jump", i.e. the closest to what you need. You can play with different values to see the effect. Note, when $x<p$ the denominator is a big number and so the fraction become close to zero. when $x>p$ the exponential in the denominator is negligible and therefore the fraction is almost the nominator.

Thus, this way you solve a problem without singularities or jumps. However, you solve a little bit different problem, which is asymptotically (as $M\to\infty$) tends to the original one.

• What condition on an interface could we pick for example? That is the point where I am facing difficulties...  I have to use finite difference methods to find the approximation of the solution. Nov 15 '15 at 22:52
• edited the post a little. the choice of interface condition depends on the problem. all the methods I've mentioned are based on FD. Nov 15 '15 at 23:06
• Could you explain me why we can take the conditions you mentioned above for the points on the interface? Nov 15 '15 at 23:12
• @evinda - all codes in one big pastebin! It seems you are not calculating the error correctly; you want $\epsilon=\sqrt{\sum_{i=1}^{N^{2}} \left(u_{e,i}-u_{a,i}\right)^2}$. In matlab you calculate that using udif = uexact-uapprox; error=norm(udif(:));. That properly calculates the 2-norm of all data points. norm(uexact-uapprox,Inf) calculates the largest singular value which is not the same. You can verify this: A = magic(3); n = sqrt(sum(sum(A.^2))), n == norm(A(:)), n == norm(A,Inf), n==norm(A) -> the first conditional returns true, the second and third conditional return false Nov 20 '15 at 11:59
• I think that we finished with the original question and now discussing another question(debugging your code), which should be posted separately, with a link to this one if you need. Otherwise we do endless comments which are less convenient. Nov 22 '15 at 10:28