NDSolve with a 0/0 error

Question

I am trying to solve this set of coupled first order differential equations.  The $phi'$ term gives a 0/0 error for when $epsilon=0$. When I enter:
NDSolve[{ϕ'[ϵ] == -1/ϵ^2 ((1 + (8 c ϕ[ϵ])/αm[M]) (Mhat[ϵ] + 
   c αm[M]^3 ϵ^3 ghat[ϕ[ϵ]/αm[M]]))/(
1 - ((16 c)/αm[M]) (Mhat[ϵ]/ϵ)), 
  Mhat'[ϵ] == αm[M]^3 ϵ^2 (((ϕ[ϵ]/αm[M])^2 - 1)^(3/2) + 3 c fhat[ϕ[ϵ]/αm[M]]), 
 Mhat[0] == 0, ϕ[0] == 1, ϕ'[0] == 0}, {Mhat, ϕ}, {ϵ, 0, 1}]

Where:
ghat[y_] := 3 ArcSinh[Sqrt[y^2 - 1]] + y (2 y^2 - 5) Sqrt[y^2 - 1]
fhat[y_] := y (2 y^2 - 1) Sqrt[y^2 - 1] - ArcSinh[Sqrt[y^2 - 1]]

and $alpha m[M]$ is just a constant.
I get
Power::infy: Infinite expression 1/0. encountered.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
Power::infy: Infinite expression 1/0.^2 encountered.
NDSolve::ndnum: Encountered non-numerical value for a derivative at [Epsilon] == 0.`.

I have tried adding the condition $phi'[0]==0$, but then I get:
NDSolve::icordinit: The initial values for all the dependent variables are not explicitly specified. 
NDSolve will attempt to find consistent initial conditions for all the variables.
NDSolve::ndnco: The number of constraints (3) (initial conditions) is not equal to the total differential order of the system plus the number of discrete variables (2).

I have also tried replacing the $phi'$ expression with a piecewise function that is the same as the above expression for $phi'$ except that it is set to zero at $epsilon=0$.  This gives me a solution that is only valid at $epsilon=0$.
I know this can be solved numerically, I just can't figure out how to get mathematica to do it.

bbgodfrey · Answer

This system of equations can be solved as follows.  For convenience, define
eq1 = ϕ'[ϵ] == -1/ϵ^2 ((1 + (8 c ϕ[ϵ])/αm[M]) (Mhat[ϵ] + c αm[M]^3 ϵ^3 ghat[ϕ[ϵ]
    /αm[M]]))/(1 - ((16 c)/αm[M]) (Mhat[ϵ]/ϵ));
eq2 = Mhat'[ϵ] == αm[M]^3 ϵ^2 (((ϕ[ϵ]/αm[M])^2 - 1)^(3/2) + 3 c fhat[ϕ[ϵ]/αm[M]]);

Because the system of equations is singular at ϵ = 0, expand the right side of eq2 there.
Normal@Series[Last[eq2], {ϵ, 0, 2}] /. ϕ[0] -> 1
(* ϵ^2 αm[M]^3 (((1 - αm[M]^2)/αm[M]^2)^(3/2) + 3 c (-ArcSinh[Sqrt[-1 + 1/αm[M]^2]] + 
   ((-1 + 2/αm[M]^2) Sqrt[(1 - αm[M]^2)/αm[M]^2])/αm[M])) *)

It follows that MHat varies as ϵ^3 near ϵ = 0, so make the substitution, Mhat[ϵ] -> ϵ^3 Nhat[ϵ] and expand the right side of eq1.
Normal@Series[Last[eq1 /. Mhat[ϵ] -> ϵ^3 Nhat[ϵ]], {ϵ, 0, 1}] /. ϕ[0] -> 1
(* -ϵ (1 + (8 c)/αm[M]) (Nhat[0] + c αm[M]^3 (3 ArcSinh[Sqrt[-1 + 1/αm[M]^2]] + 
   ((-5 + 2/αm[M]^2) Sqrt[(1 - αm[M]^2)/αm[M]^2])/αm[M])) *)

which now is well behaved and also consistent with the OP's desire that ϕ'[0] == 0.  So, make the Nhat substitution into the two equations.
eq1r = Simplify[eq1 /. Mhat[ϵ] -> ϵ^3 Nhat[ϵ]]
(* ϕ'[ϵ] == 1/(16 c ϵ^2 Nhat[ϵ] - αm[M]) ϵ (αm[M] + 8 c ϕ[ϵ]) (Nhat[ϵ] + 
   c (3 ArcSinh[Sqrt[-1 + ϕ[ϵ]^2/αm[M]^2]] αm[M]^3 - 5 αm[M]^2 ϕ[ϵ] 
   Sqrt[-1 + ϕ[ϵ]^2/αm[M]^2] + 2 ϕ[ϵ]^3 Sqrt[-1 + ϕ[ϵ]^2/αm[M]^2])) *)

eq2r = Simplify[((#/ϵ^2) & /@ eq2) /. Mhat -> Function[{ϵ}, ϵ^3 Nhat[ϵ]]]
(* 3 Nhat[ϵ] + 3 c ArcSinh[Sqrt[-1 + ϕ[ϵ]^2/αm[M]^2]] αm[M]^3 + 
   3 c ϕ[ϵ] (αm[M]^2 - 2 ϕ[ϵ]^2) Sqrt[-1 + ϕ[ϵ]^2/αm[M]^2] + 
   ϵ Nhat'[ϵ] == αm[M]^3 (-1 + ϕ[ϵ]^2/αm[M]^2)^(3/2) *)

In order that eq2r be well behaved at ϵ = 0, the boundary condition Nhat[0] must be chosen such that eq2r reduces to 0 == 0 there.
Nhat0 = Nhat[0] /. Flatten@Solve[eq2r /. ϵ -> 0 /. ϕ[0] -> 1, Nhat[0]]
(* 1/3 (6 c Sqrt[-1 + 1/αm[M]^2] + Sqrt[-1 + 1/αm[M]^2] αm[M] - 
   3 c Sqrt[-1 + 1/αm[M]^2] αm[M]^2 - 3 c ArcSinh[Sqrt[-1 + 1/αm[M]^2]] αm[M]^3 - 
   Sqrt[-1 + 1/αm[M]^2] αm[M]^3) *)

With these transformations the ODEs can be solved by
ϵ0 = 10^-6;
s = NDSolveValue[{eq1r, eq2r, ϕ[ϵ0] == 1, Nhat[ϵ0] == Nhat0} /. {c -> 1, αm[M] -> 1/6}, 
    {ϕ[ϵ], Nhat[ϵ]}, {ϵ, ϵ0, 1}];
Plot[First@s, {ϵ, ϵ0, 1}, PlotRange -> {0, 1}, AxesLabel -> {ϵ, ϕ}, 
    ImageSize -> Large, LabelStyle -> {15, Bold, Black}]
Plot[ϵ^3 Last@s, {ϵ, ϵ0, 1}, AxesLabel -> {ϵ, Mhat}, ImageSize -> Large, 
    LabelStyle -> {15, Bold, Black}]

In the absence of definitions in the Question for the two constants, I chose c = 1 arbitrarily, and αm[M] slightly smaller than ϕ[1], so that -1 + ϕ[ϵ]^2/αm[M]^2 > 0.

NDSolve with a 0/0 error

One Answer

Add your own answers!

Ask a Question