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

$\endgroup$
  • $\begingroup$ what is Epsilon? does it have value? $\endgroup$ – Algohi Sep 14 '14 at 23:38
  • $\begingroup$ @Algohi Note the last line, Block[{Epsilon = $MachineEpsilon}... $\endgroup$ – Michael E2 Sep 14 '14 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 '14 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 '14 at 4:05
  • $\begingroup$ @xzczd How can I look for a proper initial guess? Thanks $\endgroup$ – marcoac14 Sep 15 '14 at 16:11
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...

$\endgroup$
  • $\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 '14 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 '14 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 '14 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 '14 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 '14 at 1:35

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.