1
$\begingroup$

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.

Solution doesn't meet boundary condition

Improve this question
$\endgroup$
  • 2
    $\begingroup$ Are you suggesting that the coefficients $A,B,C,D$ depend on $f$? Because the way you write the ODE would imply that it's linear. $\endgroup$ – Wolfgang Bangerth Aug 24 '18 at 3:16
  • $\begingroup$ @WolfgangBangerth, I was being thick, and you're absolutely right -- this equation is linear after all. (I should have said, "A, B, C, D are nonlinear functions of x," which obviously isn't the same thing.) $\endgroup$ – emprice Aug 24 '18 at 18:41
  • $\begingroup$ Then the equation should be well-posed unless $D(x)=0$ at some points $x$. If $D$ never changes sign and is, ideally, continuous, then there should not be a problem. $\endgroup$ – Wolfgang Bangerth Aug 24 '18 at 20:17
  • $\begingroup$ @WolfgangBangerth, suppose $D(x)$ can change signs (I'm honestly not sure either way, but I think it could), but is continuous. What kinds of problems could I expect to encounter? $\endgroup$ – emprice Sep 1 '18 at 23:48
  • $\begingroup$ Then you're in trouble -- you have an ODE that is posed on two parts of the domain separately but that are inconnected or connected only in some sense that requires additional conditions. $\endgroup$ – Wolfgang Bangerth Sep 3 '18 at 13:42
1
$\begingroup$

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).

Improve this answer
$\endgroup$
  • $\begingroup$ Sounds like you talk about very specific types of boundary conditions, as well as very specific ways of incorporating them. $\endgroup$ – Christian Waluga Aug 24 '18 at 6:10
  • $\begingroup$ I just wanted to say that for every finite element or finite difference type methods, as long as you assemble their stiffness matrix (note: it's just a name for this type of linear matrix equations: $A \mathbf{x} = b$, where $A$ is called stiffness or mass matrix, $\mathbf{x}$ is the unknown vector, and $b$ is your response vector which could be assembled regarding to your boundary conditions and main ODE), you should get your imposed boundary conditions correct. Because the boundary conditions are the things that you imposed on the system and it's not related to the ODE itself. $\endgroup$ – GGG Aug 24 '18 at 13:03
  • $\begingroup$ Of course, boundary conditions should make sense if you are trying to solve a physical problem, but from pure mathematical point of view, as long as you put your intended values for boundary conditions in the $A$ and $b$ you will get it correctly at vector $\mathbf{x}$. Otherwise you did not solve that linear equation correctly, which is a basis of finite element or finite difference type methods. $\endgroup$ – GGG Aug 24 '18 at 13:05
  • 1
    $\begingroup$ If you impose boundary conditions weakly rather than strongly, for instance, I would not be so sure that it is “impossible” for the formulation to violate the boundary conditions. $\endgroup$ – Christian Waluga Aug 25 '18 at 10:31
  • 2
    $\begingroup$ See for instance Nitsche-type methods, DG or some mixed finite element methods. Of course, all methods, if correctly implemented, will have the boundary conditions built-in in some discrete fashion, but it is not always as easy to spot as in textbook examples where you end up with a bunch of equations basically saying ui=g at the boundary nodes. And, if you mess up for instance stabilization parameters or violate inf-sup conditions, etc., complete bogus may be produced at the boundaries. $\endgroup$ – Christian Waluga Aug 26 '18 at 15:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.