I've got a set of equations of motion here (describing the interior of a charged black hole, if you're interested). It's just four first order ODEs, so Mathematica should be able to solve this easily but I'm not really sure how to interpret the errors it's giving me.

Here is my code:

Q = .0025;
M = 1;
bi = .127;
xi = Q/(4*3.1416);

deqns = {eps1'[T] == 32*3.1416/bi (eps3[T]*a3[T]/(eps1[T]) + a1[T]), 
         eps3'[T] == 16*3.1416 (eps3[T]*a1[T]/eps1[T]),
         a1'[T] == 32*3.1416/bi*(eps3[T]*a3[T]*a1[T]/(eps1[T])^2),
         a3'[T] ==  16*3.1416 (a3[T]*a1[T]/(eps1[T]*bi) -  3.1416*bi*(xi^2 + 4)/eps3[T]^2)};

T0 = M + (M^2 - Q^2)^.5 - .0001;
Tf = M - (M^2 - Q^2)^.5 + .0001;

B0 = 2 M/T0 - 1 - Q^2/T0^2;
Bdot0 = 2 Q^2/T0^3 - 2 M/T0^2;
N0 = T0^2/(Abs[T0^2*B0])^.5;

ics = {eps3[T0] == T0^2,
       eps1[T0] == (T0^2*B0)^.5,
       a1[T0] == -bi/N0,
       a3[T0] == -bi*Bdot0/(2*N0*(B0)^.5)};

sol = NDSolve[{deqns, ics}, {eps1, eps3, a1, a3}, {T, Tf, T0}]

To which Mathematica replies with:

NDSolve::nlnum: The function value {112.87,5.21436,-3554.23,6.38373 +16 pi} is not a list of numbers with dimensions {4} at {T, a1[T], a3[T], eps1[T], eps3[T], (a1^[Prime])[T], (a3^[Prime])[T], (eps1^[Prime])[T], 16 pi} = {1.9999, -.000449047, 0.0158774, 0.0141418, 3.99959, 0., 0., 0., 0.}.

NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. >>

These error messages are disturbing to me for two reasons:

  1. I'm not sure why there's a 16 pi in there. Whenever I put in π, I put it in as 3.1516, not as a symbol like what appears in that error message.

  2. It's asking that I "Try giving initial conditions for both values and derivatives of the functions" but these are first order differential equations! I should not need to give initial values of the derivative, only the function.

Can anyone spot what went wrong? Am I just using NDSolve wrong?

  • $\begingroup$ It's working fine for me. You might have a pre-defined value somewhere. Restart your Kernel. I get solutions (except a singularity at T = 1.9998948863867934 that you might be able to fix with different methods) $\endgroup$
    – MathX
    Apr 7, 2016 at 0:37
  • $\begingroup$ Oh, awesome! I just copy pasted into a new worksheet and it works. I'm new to mathematica, is there an easier way to restart my Kernel? The singularity at T=1.9998948863867934 is fine, that's the event horizon of the black hole $\endgroup$
    – Mason
    Apr 7, 2016 at 0:50
  • $\begingroup$ @Mason You can use the menu: "Evaluation" -> "Quit Kernel" -> "Local" will quit the kernel. You can also use the Quit[] command. $\endgroup$
    – MarcoB
    Apr 7, 2016 at 0:54
  • $\begingroup$ @Mason Is it fine? Since NDSolve is integrating from the outside in, the domain of integration is essentially nil. $\endgroup$
    – bbgodfrey
    Apr 7, 2016 at 0:55
  • $\begingroup$ Well there are other ways, I know of Clear and Remove. Maybe read the "Lingering Definitions: when calculations go bad" part of this link: mathematica.stackexchange.com/questions/18393/… (btw this link is awesome, read it!) $\endgroup$
    – MathX
    Apr 7, 2016 at 0:56