Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to solve very simple system of coupled differential equations using the NDSolve function, but I am getting "Infinite expression 1/0.^(1/4) encountered."

Any idea how to solve the problem? I was looking for solutions where people suggested to shift the "Initial boundary point" or by using the Method -> {"Shooting", "StartingInitialConditions" ->...... Could not manage to resolve it.

Thanks in advance.

The code is as follows:

alphaphi = 2*10^-13;
MX = 1.15*10^-12;
Mphi = 10^13;
alphaX = 10^-2;
gstar = 200;
HubbleI = Mphi;
Mpl = 1.22*10^19;

phiI = (3/(8*Pi))*Mpl^2/Mphi^2 *HubbleI^2/Mphi^2;

T[x] = (30/(gstar*Pi^2))^(1/4)*Mphi*R[x]^(1/4)/x;
Xeq[x] = MX^3/Mphi^3 *(2*Pi)^(-3/2)*x^3*(T[x]/MX)^(3/2)*
Exp[-MX/T[x]];
c1 = Sqrt[3/(8*Pi)]*Mpl*alphaphi/Mphi;
c2 = c1*Mphi*alphaX/(MX*alphaphi);
c3 = c2*Mphi/MX;

sol = NDSolve[{D[phi[x], x] == -c1*x*phi[x]/Sqrt[phi[x]*x + R[x]], 
D[R[x], x] == 
 c1*x^2*phi[x]/Sqrt[phi[x]*x + R[x]] + 
  c2*x^-1*(X[x]^2 - Xeq[x]^2)/Sqrt[phi[x]*x + R[x]], 
D[X[x], x] == -c3*x^-2*(X[x]^2 - Xeq[x]^2)/Sqrt[phi[x]*x + R[x]], 
phi[1] == phiI, R[1] == 0, X[1] == 0}, {phi, R, X}, {x, 1, 10^5}];

Plot[Evaluate[{phi[x], R[x], X[x]} /. sol], {x, 1, 10^5}, 
PlotRange -> All]
share|improve this question

1 Answer 1

I suppose you know that you're dividing by R[x]^(1/4), since it's your equations. And therefore I assume the question is how to integrate the system. Since the division occurs inside an exponential and the argument of the exponential approaches negative infinity at R[1] == 0, the term stays bounded and approaches zero. One solution is to rewrite the Xeq function so that it is zero when R is zero. We also have to rewrite the T function.

Clear[Xeq, T];
T[x_, r_] := (30/(gstar*Pi^2))^(1/4)*Mphi*r^(1/4)/x;
Xeq[_, r_ /; r == 0] = 0.;
Xeq[x_?NumericQ, r_] := 
  MX^3/Mphi^3*(2*Pi)^(-3/2)*x^3*(T[x, r]/MX)^(3/2)*Exp[-MX/T[x, r]];

To get past x == 1, we need to set the option Method -> "StiffnessSwitching".

alphaphi = 2*10^-13;
MX = 1.15*10^-12;
Mphi = 10^13;
alphaX = 10^-2;
gstar = 200;
HubbleI = Mphi;
Mpl = 1.22*10^19;

phiI = (3/(8*Pi))*Mpl^2/Mphi^2*HubbleI^2/Mphi^2;

Clear[Xeq, T];
T[x_, r_] := (30/(gstar*Pi^2))^(1/4)*Mphi*r^(1/4)/x;
Xeq[_, r_ /; r == 0] = 0.;
Xeq[x_?NumericQ, r_] := 
  MX^3/Mphi^3*(2*Pi)^(-3/2)*x^3*(T[x, r]/MX)^(3/2)*Exp[-MX/T[x, r]];

c1 = Sqrt[3/(8*Pi)]*Mpl*alphaphi/Mphi;
c2 = c1*Mphi*alphaX/(MX*alphaphi);
c3 = c2*Mphi/MX;

eqns = {D[phi[x], x] == -c1*x*phi[x]/Sqrt[phi[x]*x + R[x]], 
   D[R[x], x] == 
    c1*x^2*phi[x]/Sqrt[phi[x]*x + R[x]] + 
     c2*x^-1*(X[x]^2 - Xeq[x, R[x]]^2)/Sqrt[phi[x]*x + R[x]], 
   D[X[x], x] == -c3*
     x^-2*(X[x]^2 - Xeq[x, R[x]]^2)/Sqrt[phi[x]*x + R[x]], 
   phi[1] == phiI, R[1] == 0, X[1] == 0};

sol = NDSolve[eqns, {phi, R, X}, {x, 1, 10^5}, 
   Method -> "StiffnessSwitching"];

Plots, scaled to similar ranges:

Plot[
 {phi[x]/phi[10^5], R[x]/R[10^5], X[x]/X[0.99 * 10^5]} /. sol // Evaluate,
 {x, 1, 10^5}]

Mathematica graphics

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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