I have a problem to resolve with the Finite Difference method in $[a,b]$: $$-\frac{d}{dx}(\alpha(x)\frac{du}{dx})= g(x),$$ with $\alpha(x) \in L^{\infty}$ continuous in $]a,c[$ and $]c,b[$ and discontinuous at $c$, $\alpha_{*}< \alpha(x)<\alpha_{**}$, $g \in L^2$, and $u \in C^4([a,b])$.

I want to resolve it with the FDM, I obtain a scheme on $]a,c[$ and $]c,b[$ but I don't know how to deal with the discontinuity at $c$. I take for my discretization $hc_{-}=\frac{c-a}{M+1}$ and $hc_{+}=\frac{b-c}{N+1}$ with $x_{i}^{-} = a+ihc_{-}$ and $x_{j}^{+}=c+jhc_{+}$. So I have for $1< i < M-1$ :$-\frac{\alpha_{i+1}(u_{i+1}-u_{i})}{hc_{-}^2}$ $\frac{-\alpha_{i-1}(u_{i}-u_{i-1})}{hc_{-}^2}$ for $]a,c[$ and for $M+1<i<N+M$: $-\frac{\alpha_{i+1}(u_{i+1}-u_{i})}{hc_{+}^2}$ $\frac{-\alpha_{i-1}(u_{i}-u_{i-1})}{hc_{+}^2}$ for $]c,b[$.

But i have no idea how to do for $x_M$, the point of discontinuity of the function $a$... Can you help me please ?

  • $\begingroup$ I don't know enough details to constitute a complete answer, but the immersed interface method is one approach to such problems. C.f. Leveque & Li 1994 doi.org/10.1137/0731054 or Li 2003 doi.org/10.11650/twjm/1500407515 $\endgroup$
    – whpowell96
    Commented Nov 25, 2023 at 17:47
  • $\begingroup$ That's all I have as guidance for this exercise, $\alpha_*$ and $\alpha_{**} \in R_{+}^+$ $\endgroup$
    – Kaneki Ken
    Commented Nov 25, 2023 at 19:38
  • 3
    $\begingroup$ I will note that you will have to have a very specific $g$ so that you end up with $u\in C^4$ despite the fact that the coefficient has a discontinuity. $\endgroup$ Commented Nov 26, 2023 at 3:35
  • $\begingroup$ I take it in $C^4$ for Taylor expression with "hc-" and "hc+" $\endgroup$
    – Kaneki Ken
    Commented Nov 26, 2023 at 10:28
  • 1
    $\begingroup$ @nicoguaro i pose $\varphi(x) = \alpha(x)\frac{du}{dx}$ and $\varphi' \in L^2$ $\endgroup$
    – Kaneki Ken
    Commented Nov 27, 2023 at 10:20

1 Answer 1


Since this is apparently a homework problem, let's just illustrate the idea on a simple small example. Let's take the domain [0,1], with the discontinuity at $x=0.5$, and assume $\alpha$=1 to the left of $x=0.5$, and $\alpha=2$ to the right; also assume constant $g$. In addition, assume Dirichlet boundary conditions $u(0)=u_L$ and $u(1)=u_R$

Our ODE becomes $ u_{xx} = g $ in the left half-domain, and $ u_{xx} = g/2 $ in the right half-domain. Next, to derive the matching condition at the discontinuity, integrate the ODE over a narrow segment covering the discontinuity,

$ \int_{0.5-\epsilon}^{0.5+\epsilon} \left[ \frac{d}{dx} ( \alpha u') = g \right] dx, $

which leads to

$ ( \alpha u')_{0.5+\epsilon} - ( \alpha u')_{0.5-\epsilon} = 2 g \epsilon, $

and in the limit $\epsilon \rightarrow 0$ for our choice of $\alpha(x)$ it becomes

$u'_{x-} = 2 u'_{x+}$,

using the one-sided derivative notation; this matching condition has the obvious meaning of flux conservation.

Now, let's do FD discretization of the ODE using just five uniformly distributed grid points and using simplest FD approximations; $u_0$ and $u_4$ correspond to the boundaries, and $u_2$ corresponds to $x=0.5$. This FD formulation amounts to solving the following system of five linear equations:

$u_0 = u_L$

$u_0 - 2 u_1 + u_2 = h^2 g $

$2(u_3-u_2) = u_2-u_1$

$u_2 - 2 u_3 + u_4 = h^2 g/2 $

$u_4 = u_R$

Here $h$ is the grid spacing, the third equation is the matching condition at the discontinuity.

It is easy to verify that this linear system is solvable:

import numpy as np
A = np.array([[1,0,0,0,0], [1,-2,1,0,0], [0,1,-3,2,0], [0,0,1,-2,1],[0,0,0,0,1]])
  • $\begingroup$ But we still don't have the form of the approximation for the discontinuity, what do you have for f(c) ? for example, If $c = x_k$, $h^2g_k = ... ?$, I think for the example : 2u_1 - 3u_2 - 2u_3 = f_k ??? $\endgroup$
    – Kaneki Ken
    Commented Nov 28, 2023 at 21:49
  • $\begingroup$ Well, if the RHS source function g does not contain a delta-function at the discontinuity, we have the flux conservation across the discontinuity, $(\alpha u')_{-} = (\alpha u')_{+}$. That's the equation describing the discontinuity. $\endgroup$ Commented Nov 29, 2023 at 0:12
  • $\begingroup$ so what I wrote for the nodes of the discontinuity is good?, and for the generalization I think I have to use $\alpha^-$ and $\alpha^+$ at each side and do the subtraction no? $\endgroup$
    – Kaneki Ken
    Commented Nov 29, 2023 at 7:16
  • 1
    $\begingroup$ @KanekiKen For the gird node corresponding to the discontinuity, you'll need to use that special equation expressing the flux conservation $(\alpha u')_{-} = (\alpha u')_{+}$, that's it. $\endgroup$ Commented Nov 30, 2023 at 17:17

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.