4

Hey guys I need someone to give me a hand on this. I don't know if this is too complicated or it's just the lack of knowledge I have on Mathematica.

I'm trying to solve the following equation but NDSolve triggers several errors.

data = {Q -> 1, R -> 1}

(* Out[36]= {Q -> 1, R -> 1} *)

eqs = {w''[x] w'[x] + 3 Q w[x]^2 - 2 R w[x] w'[x] == 0, w[Epsilon] == 1, w[1] == 0.8} /. data;

s = Block[{Epsilon = $MachineEpsilon}, NDSolve[eqs, w, {x, Epsilon, 1}]]

Errors triggered

Power::infy: "Infinite expression 1/0. encountered."
Infinity::indet: "Indeterminate expression 0.\ ComplexInfinity encountered."
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 2.220446049250313`*^-16.

Thanks in advance

Michael E2
  • 235,386
  • 17
  • 334
  • 747
Marco
  • 41
  • 2

1 Answers1

3

I think the problem is not with Mathematica or your understanding of it but with the mathematics of that equation and it isn't even a singularity which would show for Epsilon->0. The problem is that your equations have a problem when the derivative w'[x] becomes 0, and that does happen already for much higher values of Epsilon (I have renamed that to eps in the following):

data = {Q -> 1, R -> 1}
eqs[eps_] := {
      w''[x] w'[x] + 3 Q w[x]^2 - 2 R w[x] w'[x] == 0, 
      w[eps] == 1, w[1] == 0.8
  } /. data;
eps = 0.75567;
s = NDSolveValue[eqs[eps], w, {x, eps, 1}];
Plot[{s[x], s'[x]}, {x, eps, 1},Frame->True]

enter image description here

note that the critical point isn't at the lowest x-values but at the endpoint x=1, for roughly the given eps-value the derivative effectively gets zero and that seems to be the source of your problem. Solving the differential equations for w''[x] is another way to see the source of the problem. Here is an Animation that illustrates what happens when you decrease the x value where the initial condition is set:

Animate[
  s = NDSolveValue[eqs[eps], w, {x, eps, 1}];
  Plot[{s[x], s'[x]}, {x, eps, 1}, Frame -> True, 
  PlotRange -> {{0.7, 1}, {-4, 2}}],
  {eps, 0.75567, 0.99}, AnimationDirection -> Backward
]

EDIT Actually you will find that older versions of Mathematica won't return a result even for values for eps larger than the one I have given. This does indicate that the equations are problematic and as has been mentioned by xzczd that is of course an indication that the results that Version 10 gives might also not be trusted blindly. If you know that there are sigularities in a differentital equations any way to solve it numerically needs some extra effort to justify that the result is correct.

EDIT In a comment xzczd has mentioned that the original code didn't return the same result reliably when playing around with eps. I have also seen that behavior and updated the code in a way that I think will not show that problem. It is a good example why it is always a good idea to make dependendcies of your equations/formulae/functions explicit...

Albert Retey
  • 23,585
  • 60
  • 104
  • Thanks for helping. What's the output of this code? I tried to evaluate it but it still triggers the same errors. – Marco Sep 16 '14 at 01:40
  • @Marco I think Albert is in v10. Just tested the code on the wolfram cloud and it does work. – xzczd Sep 16 '14 at 03:28
  • @AlbertRetey I think the v10 result isn't quite reliable. The same code can produce different result. For example, make the code fail by setting eps = 0.35 etc. and try eps = 0.75567 again, this time the code will fail. In some conditions I can even make it success but give different plots! – xzczd Sep 16 '14 at 03:42
  • 1
    @xzczd: I absolutely aggree that one shouldn't trust results of numeric algorithms in such situations (known potential singularities) blindly, and I have admittedly not made any tests about the quality of what NDSolve returns (and yes, only v10 will return a result at all). In this case I have also seen what you have seen, but I think that is because the current formulation is suboptimal in how it handles the eps parameter. I will update to show a more reliable way... – Albert Retey Sep 16 '14 at 14:29
  • @AlbertRetey I'm not sure if I got your point. You said that the critical point is at x=1 and that's where my problem lies. But I don't get why it only triggers errors when eps decreases. My professor suggested me to add a transient term so I can emulate a pseudo-transient. However I need to estimate w[x,0] so it converges faster. What do you think about this solution? Yes, I can estimate w[x,0] using a simplified model. – Marco Sep 17 '14 at 01:35
  • @xzczd I thought I had tested on version 10, but I hadn't. I just made it work, but I'm facing the issue you guys have faced. Thanks for helping – Marco Sep 17 '14 at 01:41
  • @Marco: when you change eps you solve a different problem: if eps is > 0.75567 the solution to your problem will have w'[x]<0 for 0.75567 < x < 1 and there is no problem with the singularity. If eps gets about 0.75567, then the solution (at least the one Mathematica finds/tries with the shooting methods) will have w'[1]==0 and hits the singularity. It could be that there are other solutions which the shooting method as implemented in Mathematica just doesn't find (or you have to give it an initial guess). I don't think this will work here... – Albert Retey Sep 17 '14 at 08:57
  • ... so you probably are best off to try what your professor suggested (another wellknown method to approach such problems). I would also think that it might be worth to rethink what problem that equation originated from and probably understand whether that singularity has a meaning in that context. Maybe you then can reformulate the problem in a way that is less problematic for numeric solvers... – Albert Retey Sep 17 '14 at 09:01