# My toy Laplace equation solver using finite-difference is unstable and I'm not sure why

I am trying to solve the variable-coefficient Laplace equation $$\partial(\epsilon\partial u) = (\partial\epsilon)(\partial u) + \epsilon\partial^2 u = 0$$using a finite difference scheme: $$\left(\frac{\epsilon(x+h) - \epsilon(x-h)}{2h}\right)\left(\frac{u(x+h)-u(x-h)}{2h}\right) + \epsilon(x) \frac{u(x+h) - 2u(x) + u(x-h)}{h^2}=0$$This gives us $$u(x) = \frac{1}{2}\left(\frac{\epsilon(x+h) - \epsilon(x-h)}{2}\right)\left(\frac{u(x+h)-u(x-h)}{2}\right) + \frac{\epsilon(x)}{2}\left[u(x+h) + u(x-h)\right]$$Here is my python script that solves the equation on the domain [0,1], and sets the coefficient to 1 everywhere except on the interval [0.4, 0.6] where it is set to 2. Instead of converging, it produces nans. If I set the coefficient to the same everywhere, it works fine.

import numpy as np
import matplotlib.pyplot as plt

num_points = 100
domain = np.linspace(0, 1, num_points)

# Coefficient is 1 everywhere except in the interval [0.4, 0.6]
eps = np.ones(num_points)
for i, x in enumerate(domain):
if 0.4 <= x <= 0.6:
eps[i] = 2.0

# Solution
u = np.zeros(num_points)

#  Boundary conditions: u(0) = 1, u(1) = 0
u[0] = 1.0

while True:
u_old = u.copy()
for i in range(1, num_points - 1):
u[i] = 0.125 * (eps[i + 1] - eps[i - 1]) * (u[i + 1] - u[i - 1])
u[i] += 0.5 * eps[i] * (u[i + 1] + u[i - 1])

if np.allclose(u, u_old):
break

plt.plot(domain, u)
plt.show()


The iteration could be blowing up.

You can see this by looking at your iteration

$$|u_{i+1}|\le \frac{1}{4} \omega(\epsilon, 2h) |u_i|_{\infty} + |\epsilon|_{\infty} |u_i|_{\infty}$$

where $$\omega(\epsilon, 2h) = \sup |\epsilon(x) - \epsilon(x+2h)|$$. In your homogeneous case, $$\omega = 0$$, and $$|\epsilon| = 1$$, there is no blowing up.

For discontinuous case, $$\omega = 1$$, then

$$|u_{i+1}|\le \frac{5}{4} |u_i|$$

which does not have any control over the maximum.

To verify this, if you bring down your $$\epsilon$$ to something like $$\epsilon = 0.25$$ overall except $$\epsilon = 0.5$$ on some parts, then there should be no problem.

• Thanks. Could this issue also be resolved by evaluating the coefficient on a different grid from the solution? I.e. Is an expression like this valid? mathb.in/77138 where h' > h Commented Dec 11, 2023 at 10:35
• yes. But the magnitude of $\epsilon$ maybe the bottleneck. @DJames, it seems your iteration should divide $\epsilon$ on rhs. Commented Dec 11, 2023 at 17:06
• It may also help to apply the finite-difference scheme to the conservative form of the equation, $(\epsilon u_x)_x=0$. If $\epsilon$ is not differentiable, then applying the product rule may not be justified and your discretization is having problems due to this. Commented Dec 11, 2023 at 17:23
• @whpowell96 unfortunately that is likely my problem: In my use case, the coefficient changes abruptly, like a step function, at various interfaces. Here is my original question that motivated this follow-up physics.stackexchange.com/q/791946/136054 Commented Dec 12, 2023 at 16:10
• It will definitely help to try the conservative discretization. If you still have problems, there are specific methods for interface problems such as this. C.f., The immersed interface method. Commented Dec 12, 2023 at 16:34