Finite difference problem

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

• 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 Commented Nov 25, 2023 at 17:47
• That's all I have as guidance for this exercise, $\alpha_*$ and $\alpha_{**} \in R_{+}^+$ Commented Nov 25, 2023 at 19:38
• 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. Commented Nov 26, 2023 at 3:35
• I take it in $C^4$ for Taylor expression with "hc-" and "hc+" Commented Nov 26, 2023 at 10:28
• @nicoguaro i pose $\varphi(x) = \alpha(x)\frac{du}{dx}$ and $\varphi' \in L^2$ Commented Nov 27, 2023 at 10:20

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

$$( \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]])
np.linalg.det(A)
-6.000000000000001

• 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 ??? Commented Nov 28, 2023 at 21:49
• 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. Commented Nov 29, 2023 at 0:12
• 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? Commented Nov 29, 2023 at 7:16
• @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. Commented Nov 30, 2023 at 17:17