4
$\begingroup$

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

Improve this question
$\endgroup$
8
  • $\begingroup$ what is Epsilon? does it have value? $\endgroup$
    – Basheer Algohi
    Sep 14, 2014 at 23:38
  • $\begingroup$ @Algohi Note the last line, Block[{Epsilon = $MachineEpsilon}... $\endgroup$
    – Michael E2
    Sep 14, 2014 at 23:47
  • $\begingroup$ @Algohi I looked for NDSolve::ndnum on help and that's just a suggestion to avoid singularity at 0. MachineEpsilon is a very small number, close enough to zero, but not zero. $\endgroup$
    – Marco
    Sep 15, 2014 at 0:19
  • 1
    $\begingroup$ Seems to be a task for shooting method, then you'll have a hard time looking for a proper initial guess… $\endgroup$
    – xzczd
    Sep 15, 2014 at 4:05
  • $\begingroup$ @xzczd How can I look for a proper initial guess? Thanks $\endgroup$
    – marcoac14
    Sep 15, 2014 at 16:11

1 Answer 1

Reset to default
3
$\begingroup$

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

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thanks for helping. What's the output of this code? I tried to evaluate it but it still triggers the same errors. $\endgroup$
    – Marco
    Sep 16, 2014 at 1:40
  • $\begingroup$ @Marco I think Albert is in v10. Just tested the code on the wolfram cloud and it does work. $\endgroup$
    – xzczd
    Sep 16, 2014 at 3:28
  • $\begingroup$ @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! $\endgroup$
    – xzczd
    Sep 16, 2014 at 3:42
  • 1
    $\begingroup$ @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... $\endgroup$
    – Albert Retey
    Sep 16, 2014 at 14:29
  • $\begingroup$ @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. $\endgroup$
    – Marco
    Sep 17, 2014 at 1:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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