0
$\begingroup$

I am trying to solve a set of four coupled equations through NDSolve. However, I am getting the error "step size is effectively zero; singularity or stiff system suspected.." I think this is a common error when there exists very large or very big numbers, but I can't see any such numbers in my system. So, I am worried what could be the reason ? The set of equations are given below :

eqn1 = M'[a] + 
Log[10]/Sqrt[1/3  (B[a] 10^(-3 a) + R [a] 10^(-4 a) )] (1/M[a]) ^2;
eqn2 = R'[a] + 10^a/M[a] M'[a] B[a];
eqn3 = B'[a] - 1/M[a] M'[a] B[a];
eqn4 =  s'[a] + 3 s[a] Log[10] + 
1/  (s[a]^((1/3))) B[a]/(M[a] 10^a) M'[a];
sol1 = NDSolve[{eqn1 == 0 , eqn2 == 0, eqn3 == 0, 
eqn4 == 0, M[2] == 2, R[2] == 1, 
B[2] == 0.001, s[2] == 0.1}, {M, R, B , s}, {a, 2, 8}]
$\endgroup$
1
  • 2
    $\begingroup$ Have you plotted the solution? If the solution has a singularity, it may be a feature of the system and there is nothing to do about it. (E.g. how would you get rid of the pole at $x=0$ of $y=1/x$?). The singularity may also be due to the system being unstable in the neighborhood of your IVP. Related: (39028), (232756). $\endgroup$ – Michael E2 Jan 12 at 16:45
0
$\begingroup$

Seems for all positive initial conditions (m2) for M, system has a singularity. But for negative m2, you get a full a-range solution.

ndsol[b2_, m2_, r2_, s2_] := 
  First@NDSolve[{eqn1 == 0, eqn2 == 0, eqn3 == 0, eqn4 == 0, 
B[2] == b2, M[2] == m2, R[2] == r2, s[2] == s2}, {B, M, R, s}, {a,
 2, 8}]

Plot[Evaluate[M[a] /. ndsol[.001, 2, 1, .1]], {a, 2, 
   2 + .9999*(puf[[1, 2, 1, 1, 2]] - 2)}]

Manipulate[
  Plot[Evaluate[M[a] /. (puf = ndsol[b2, m2, r2, s2] // Quiet)], {a, 2,
2 + .9999*(puf[[1, 2, 1, 1, 2]] - 2)}], {{b2, .001}, -3, 3, 
Appearance -> "Labeled"}, {{m2, -2}, -5, 5, 
Appearance -> "Labeled"}, {{r2, 1}, -5, 5, 
Appearance -> "Labeled"}, {{s2, .001}, -1, 1, 
Appearance -> "Labeled"}]

enter image description here

$\endgroup$

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.