0
$\begingroup$

Tried googling the error "step size is effectively zero; singularity or stiff system suspected". Haven't found any soln. Any help would mean a lot. Posting the full code.

w0 = -0.5; w1 = 0.05; wm = 0.01; \[CapitalOmega]m0 = 0.23; \
\[CapitalOmega]c0 = 0.78; \[Delta] = 0.01; H = 72;

h1[z_] := \[CapitalOmega]m0 ((1 + z)^(3 (1 + wm - \[Delta])) + 
      3 \[Delta] (1 + z)^(-3 (1 + w0 - w1))*
       E^(w1 (1 + z))*(-3 w1)*(1 + z)^(-3 (w0 - w1 - wm + \[Delta]) - 
          1)*(-(-3 w1)*(1 + z))^(3 (w0 - w1 - wm + \[Delta]) - 
          1) Gamma[-(3 (w0 - w1 - wm + \[Delta]) - 1) - 
         1, -(-3 w1) (1 + z)]) + \[CapitalOmega]c0 (1 + z)^(
    3 (1 + w0 - w1)) E^(3 w1 (1 + z));
s1 = NDSolve[{E^(K[z]/3) K'[z] == 3/Sqrt[h1[z]], K[.01] == .01}, 
   K[z], {z, 0, 10}];
f1[z_] := Evaluate[K[z] /. s1];

w00 = -1; 
h2[z_] := \[CapitalOmega]m0 ((1 + z)^(3 (1 + wm - \[Delta])) + 
      3 \[Delta] (1 + z)^(-3 (1 + w00 - w1))*
       E^(w1 (1 + z))*(-3 w1)*(1 + 
          z)^(-3 (w00 - w1 - wm + \[Delta]) - 
          1)*(-(-3 w1)*(1 + z))^(3 (w00 - w1 - wm + \[Delta]) - 
          1) Gamma[-(3 (w00 - w1 - wm + \[Delta]) - 1) - 
         1, -(-3 w1) (1 + z)]) + \[CapitalOmega]c0 (1 + z)^(
    3 (1 + w00 - w1)) E^(3 w1 (1 + z));
s2 = NDSolve[{E^(K[z]/3) K'[z] == 3/Sqrt[h2[z]], K[.01] == .01}, 
   K[z], {z, 0, 10}];
f2[z_] := Evaluate[K[z] /. s2];

w000 = -1.5;
h3[z_] := \[CapitalOmega]m0 ((1 + z)^(3 (1 + wm - \[Delta])) + 
      3 \[Delta] (1 + z)^(-3 (1 + w000 - w1))*
       E^(w1 (1 + z))*(-3 w1)*(1 + 
          z)^(-3 (w000 - w1 - wm + \[Delta]) - 
          1)*(-(-3 w1)*(1 + z))^(3 (w000 - w1 - wm + \[Delta]) - 
          1) Gamma[-(3 (w000 - w1 - wm + \[Delta]) - 1) - 
         1, -(-3 w1) (1 + z)]) + \[CapitalOmega]c0 (1 + z)^(
    3 (1 + w000 - w1)) E^(3 w1 (1 + z));
s3 = NDSolve[{E^(K[z]/3) K'[z] == 3/Sqrt[h3[z]], K[.01] == .01}, 
   K[z], {z, 0, 10}];
f3[z_] := Evaluate[K[z] /. s3];

Plot[{f1[z], f2[z], f3[z]}, {z, 0, 5}] ```
$\endgroup$
5
  • $\begingroup$ What do the last two lines have to do with your problem? Why are they included if they don't matter? $\endgroup$
    – Michael E2
    Nov 7, 2021 at 12:45
  • $\begingroup$ Didn't quite get your question? Whatever lines I have added all are necessary. Please try to run the code in your machine. P.S- Thank you for yesterday's help. It worked at last. $\endgroup$ Nov 7, 2021 at 12:49
  • $\begingroup$ The last two lines of code are f3[z_] := Evaluate[K[z] /. s3]; and Plot[{f1[z], f2[z], f3[z]}, {z, 0, 5}]. The problem occurs in the 3rd-to-last line with NDSolve. I don't want to Plot a failed NDSolve. It potentially wastes time, distracts from the real problem, causes lots of error messages, all of which are pointless until the NDSolve code is fixed. I sometimes just skip such questions. -- P.S. You're welcome. $\endgroup$
    – Michael E2
    Nov 7, 2021 at 13:03
  • $\begingroup$ Actually, you should just include the one NDSolve that fails, not all three. Why make it so difficult for other to find your problem? $\endgroup$
    – Michael E2
    Nov 7, 2021 at 13:05
  • $\begingroup$ Will keep that in mind from next time. In my machine I am getting problem with the 2nd and 3rd NDSolve. First one worked fine and gave the graph. $\endgroup$ Nov 7, 2021 at 13:09

2 Answers 2

1
$\begingroup$

The integrand in your second integration seems to be stiff what would require extremeliy small steps. Therefore, use the method "StiffnessSwitching" in the second integration like:

s2 = NDSolve[{E^(K[z]/3) K'[z] == 3/Sqrt[h2[z]], K[.01] == .01}, 
  K[z], {z, 0, 10}, Method -> "StiffnessSwitching"]

enter image description here

$\endgroup$
3
  • $\begingroup$ Tried adding ``` Method -> "StiffnessSwitching"```, the error got way but still I am getting one graph instead of three. $\endgroup$ Nov 7, 2021 at 12:51
  • $\begingroup$ @DarthVader Don't think that might be because of the same reason as last time? (Complex numbers.) Also you might want to define f1[z_] := Evaluate[K[z] /. First@s1] $\endgroup$
    – Michael E2
    Nov 7, 2021 at 13:10
  • $\begingroup$ It seems I got the same complex number problem like yesterday. Will try to rectify it my changing the parameters. $\endgroup$ Nov 7, 2021 at 13:49
1
$\begingroup$

First, don't use K as your own symbol. It is used internally by Mathematica.

You can determine whether it's stiffness (often fixable) or a singularity (a feature which does not need fixing, and can't be "fixed"). Plotting the solution or the norm of the solution might suggest a singularity. Trying a non-stiff solver such as an explicit RK method will test for stiffness.

s2 = NDSolve[{E^(k[z]/3) k'[z] == 3/Sqrt[h2[z]], k[.01] == .01},
 k, {z, 0, 10}, Method -> "ExplicitRungeKutta"]

(* {{k -> InterpolatingFunction[{{0., 10.}}, "<>"]}} *)

No stiffness error, so it must be a problem in the (Automatic) LSODA method.


Well, duh, there's an obvious singularity when h2[z] == 0, which occurs at the point z == 0.8... in the error. I should have checked this first, but a little slow this morning, I guess.

Plot[h2[z], {z, 0, 1}]

Since it's a pole of order 1/2, it's integrable, but I don't know of a graceful way to get NDSolve to solve it. The method "ExplicitRungeKutta" must accidentally skip over it due to round-off error. In fact, if we set WorkingPrecision -> 16, we get a divide-by-zero error.

$\endgroup$
1
  • $\begingroup$ We are allowed to change the parameters w0, w1, wm, $\delta$. Will see if changing them solves the issue. Meanwhile, if you were able to find something do kindly share. Thank you. $\endgroup$ Nov 7, 2021 at 13:54

Your Answer

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

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