Numerical methods for boundary-value ODEs with a jump condition

Question

I want to solve a non linear system of equations of a particular kind. I find it hard to formulate clearly so I directly give a simple example.

$ f''=A(f,g)\\ g''=B(f,g) $

with the boundary conditions as follow:

$f'(0)=g'(0)=0$

$f'(\infty)=g'(\infty)=0$

At $x=1$, f is continuous but there is a discontinuity for $g$:

$ f'(1^-)=f'(1^+), ~~~ f(1^-)=f(1^+) \\ g(1^-)=g_{left}, ~~~g(1^+)=g_{right} \\ $

So we can see this as two problems, one for $x<1$, and one for $x>1$, but the boundaries conditions are coupled in $x=1$. When I have an analytical solution, there is no big deal. Solve the problem for each side, then match the solution at $x=1$.

I'm not used with numerical methods and I don't even find how to call this problem. I have some ideas to tackle this but I wish to read a bit about that before.

Answer 1

These problems are usually called interface problems with jump conditions; an example would be the paper Immersed-interface finite-element methods for elliptic interface problems with nonhomogeneous jump conditions.

Typically, they come up in PDEs in fracture mechanics to model discontinuities in stress due to cracks, but they can also come up in problems like thermal modeling or electromagnetic modeling to depict something called a "contact resistance" or basically an appreciable resistance due to a thin film.

If you know the location of the discontinuity, and it does not move or change in shape, you would typically fit a discretization to the discontinuity, so that the discontinuity coincides with nodes on your mesh. If you duplicate mesh nodes along the discontinuity -- in your case, it's 1-D, so you would have two nodes at $x = 1$, so that $x = 1^{-}$ is one node, and $x = 1^{+}$ is another node -- then you can use that discretization to form numerical methods. This approach is one usually taken for PDEs with spatial discontinuities. You could also track the interface implicitly using cut-cell discretizations, with level sets, or some sort of front tracking.

If the discontinuity is in 1-D, adaptive ODE integration methods are also an option, and in that case, one viable strategy is to do what you propose and integrate from your initial condition to $x = 1$, stop the integrator, reinitialize it with new initial conditions at $x = 1$ that reflect the jump condition, and restart the integrator.