Can a second-order ODE be "inconsistent" with its boundary conditions?

Question

I am trying to solve a set of coupled, nonlinear ODEs. The only dependent variable is a 1-dimensional spatial coordinate, let's call it $x$. For now, I've managed to approximate away some of the coupling between the equations, such that the first one doesn't depend on the second one at all. I solved it with FEM without major difficulty: boundary conditions met, residuals look good, physically meaningful result.

So now I'm left with one nonlinear, second-order ODE where the nonlinear coefficient functions are completely determined by the FEM solution, something like

$$ \begin{cases} A(x) + B(x) f(x) + C(x) f'(x) + D(x) f''(x) = 0 \\ f(0) = f(1) = 0 \end{cases} $$

where $A$, $B$, $C$, and $D$ are known, at least at the resolution of the FEM mesh.

The problem is that, using the same basic code as before, the solution appears to violate the boundary conditions entirely, blowing up as $x \rightarrow 0$ instead of approaching zero (or even any finite value, it doesn't seem to matter what I choose). Thinking this was a problem with my code, I experimented with a finer mesh and with finite differences instead of finite elements, with no success. It finds the same nonsensical result every time.

Finally, I put my numerical data for the coefficients into Mathematica and attempted to solve the above BVP. I let it run for quite a while (maybe 15-30 minutes) before giving up. It never found a solution.

So finally, my question: Is it possible that this equation simply can't be solved with these boundary conditions? For example, I tried changing $f(1) = 0 \rightarrow f'(0) = 0$ and the solution was found almost instantaneously by Mathematica, but it still violated the new boundary condition; somehow, that didn't trigger an error, but it made me suspicious that something "bigger" might be going on here.

I'm looking for some mathematical intuition behind why the numerical methods might fail in these cases and how to mitigate these problems.

Update

This equation is linear (correcting here so that comments still make sense). I think that's not the root of my problem, though; I tried using a linear solver with the proper stiffness matrix and load vector to solve $-f''(x) = g(x)$, which should be EASY; here, $g(x)$ is the FEM solution to the first equation. The solution meets boundary conditions but is "jagged" (think sawtooth patterns overlaid on a smooth function). Now I think perhaps there's something wrong with substituting the FEM solution for the first equation into the second equation, but I'm not sure why.

Not giving up yet.

Update 2

I think I can safely rule out a problem with this code; I am generating the source files with SymPy, so there shouldn't be any human errors associated with changing functional forms, etc. of the equations I want to solve. As an example of what I'm seeing now, see figure below. The only intuition I have for this is that the code "tries" to meet boundary conditions (which are homogeneous Dirichlet), but the solution, for whatever reason, fundamentally wants to blow up to $-\infty$. So we end up with this "compromise" solution that straddles the two possibilities.

Is this even remotely correct? It doesn't seem to matter what I change at this point, and refining the mesh just increases the frequency of the sawtooth form.

Answer 1

It is impossible that finite element or finite difference violate their boundary condition that you imposed. It is not related to the physics of your problem and it's just an issue with your implementation. No matter your ODE is consistent with your boundary conditions or something else. The only thing matters is that you enforce your boundary conditions when you are trying to assemble your stiffness matrix and response vector in your finite element or finite difference methods and if you assemble them correctly it should recover the imposed boundary conditions. Otherwise maybe you have a bug in your matrix assembling which makes the stiffness matrix singular, which shows you some nonsense values at the boundaries. As a result, my suggestion is that check your matrix assembling function and also check if your linear solver or whatever could generate informative warnings/errors when something is wrong with your stiffness matrix or response vector. You need to get at least whatever you wanted to impose on your equation correctly and then that's an another story when you want to discuss about the physical meaning of your model, which makes sense or not.

Also your ODE is not nonlinear because $A(x)$, $B(x)$, $C(x)$, and $D(x)$ do not depend on $f$ itself and as a result it should be solved easily by simple linear finite element or finite difference solvers (e.g. scipy BVP solver).