Finite difference problem

Question

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 ?

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