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.$

enter image description here

$\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.

enter image description here

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$?

  • $\begingroup$ what physical problem do you considering? $\endgroup$ Nov 15, 2015 at 22:51
  • $\begingroup$ Is your interface moving? $\endgroup$
    – nluigi
    Nov 17, 2015 at 22:40
  • $\begingroup$ No, it's the line $y=0$. @nluigi $\endgroup$
    – Mary Star
    Nov 17, 2015 at 23:15
  • $\begingroup$ 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. $\endgroup$
    – nluigi
    Nov 18, 2015 at 14:04
  • $\begingroup$ @nluigi Could you look at my question of the last comment if you have an idea? $\endgroup$
    – Mary Star
    Nov 19, 2015 at 15:20

1 Answer 1


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.

  • $\begingroup$ 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. $\endgroup$
    – Mary Star
    Nov 15, 2015 at 22:52
  • $\begingroup$ edited the post a little. the choice of interface condition depends on the problem. all the methods I've mentioned are based on FD. $\endgroup$ Nov 15, 2015 at 23:06
  • $\begingroup$ Could you explain me why we can take the conditions you mentioned above for the points on the interface? $\endgroup$
    – Mary Star
    Nov 15, 2015 at 23:12
  • 1
    $\begingroup$ @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 $\endgroup$
    – nluigi
    Nov 20, 2015 at 11:59
  • 1
    $\begingroup$ 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. $\endgroup$ Nov 22, 2015 at 10:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.