2

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?

MarcoB
  • 67,153
  • 18
  • 91
  • 189
Mason
  • 560
  • 2
  • 10
  • 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) – MathX Apr 07 '16 at 00:37
  • 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

     – Mason Apr 07 '16 at 00:50
  • @Mason You can use the menu: "Evaluation" -> "Quit Kernel" -> "Local" will quit the kernel. You can also use the Quit[] command. – MarcoB Apr 07 '16 at 00:54
  • @Mason Is it fine? Since NDSolve is integrating from the outside in, the domain of integration is essentially nil. – bbgodfrey Apr 07 '16 at 00:55
  • Well there are other ways, I know of Clear and Remove. Maybe read the "Lingering Definitions: when calculations go bad" part of this link: http://mathematica.stackexchange.com/questions/18393/what-are-the-most-common-pitfalls-awaiting-new-users (btw this link is awesome, read it!) – MathX Apr 07 '16 at 00:56
  • @bbgodfrey Aw crap, you're right. This isn't behaving afterall. So, I only want it to go from T=Tf to T=T0, though my 'initial' conditions are for T=T0.

    Is that a problem in mathematica?

     – Mason Apr 07 '16 at 01:19
  • I've tried plotting eps_3[T] using the command Plot[eps3[T] /. sol, {T, Tf, T0}] and sometime is definitely wrong. eps3 should go as T^2. – Mason Apr 07 '16 at 01:21
  • @Mason This is a mathematics issue, not a Mathematica issue. You are solving an initial value problem, beginning at T0. To circumvent the apparent singularity at the event horizon, I recommend that you change coordinates to remove it, which I know is possible from my work on the Kerr metric 45 years ago. – bbgodfrey Apr 07 '16 at 01:29
  • @Mason I should add that any results you get from Plot[eps3[T] /. sol, {T, Tf, T0}] are meaningless, because sol is valid only very near T = 2. – bbgodfrey Apr 07 '16 at 01:36
  • @bbgodfrey I think you're right, something must have gone wrong in my Poisson Brackets. The coordinate transformation is not necessary since I only care about what's going on inside the horizon. I had hoped that what was wrong was just that Mathematica was evolving the solution the wrong way and hitting the coordinate singularity but that seems to not be the case. – Mason Apr 07 '16 at 01:48

0 Answers0