1
$\begingroup$

I am having some problems with handling an event during NDSolve.

The event itself is simple, I am looking for when one of my coordinates, defined in

coords = {t, r, θ, ϕ};

has a value less than some constant, which I call rplus. That is, when r

I have had a look at the reference pages for EventLocator and managed to apply it to simpler problems than mine.

For example, here I solve the 2D Keplerian motion of a planet and just to test out the EventLocator function, I make it stop the integration when the y-coordinate goes negative.

x0 = -1; y0 = 0; x0p = 0; y0p = Pi; tf = 0.5;
Aorbit = NDSolve[{x''[t] + (4 Pi^2) x[t]/(x[t]^2 + y[t]^2)^(3/2) == 
 0,
                  y''[t] + (4 Pi^2) y[t]/(x[t]^2 + y[t]^2)^(3/2) == 0,
                  x[0] == x0, y[0] == y0,
                  x'[0] == x0p, y'[0] == y0p},
                  {x, y}, {t, 0, tf},
                  Method -> {"EventLocator",
                             "Event" :> y[t] < 0,
                             "EventAction" :> Throw[tend = t, "StopIntegration"]}];

My problem involves 4 equations of this sort and sadly, I will have a hard time listing the full equations here (They are the geodesic equations for a particle moving in the Kerr spacetime).

The relevant part of the code is presented here:

deq = Table[coords[[i]]''[τ] == (geodesic[[i]] /.
      Join[Table[coords[[i]]' -> coords[[i]]'[τ], {i, 1, n}],
      Table[coords[[i]] -> coords[[i]][τ], {i, 1, n}]]), {i, 1, n}];

deq = Join[deq, Table[coords[[i]]'[0] == ivs[[i]], {i, 1, n}],
      Table[coords[[i]][0] == ics[[i]], {i, 1, n}]];

soln = NDSolve[deq, coords,{τ, 0, maxτi},
       Method -> {"EventLocator",
                  "Event" :> Re[coords[[2]][τ]] <= Re[rplus],
                  "EventAction" :> Throw[τend = τ, "StopIntegration"]}

When I run something such as (uinvar is a variable used inside computeSoln and computeSoln is a function made up of the commands above)

uinvar = -1; M = 1; a=0.9;
soln = computeSoln[10000, {0, 0, 0.0000006}, {0, 10000, π/2, Pi/4}];

I get an error message as such:

NDSolve::evre: The value of the event function at τ = 0.00010263637541427471` was not a real number. The event will be ignored in steps where it does not evaluate to real numbers at both ends

So, it seems that coords[2] is behaving badly. However, as you can see from my code, my initial guess was that coords[2] (that is, the r coordinate) has an imaginary part. And sure enough:

ir = soln[[1]][[2]][[2]]; ir[0.00001]

returns

5. + 0.*I

However, adding Re[] around both the LHS and RHS of the event didn't help.

After some searching, I found this thread where one of the answers states that one of the sides in the event is a list and not a number. However, I don't think this is the case for my problem.

So I'm stumped. I am almost certain this is some jumble due to me using lists for my coordinates and such but I can't seem to crack it. If it would help, I can try and get some minimalistic version of the whole problem such that other people could run it, please let me know if this would help.

Edit: I includ here a simple version of what I am doing, making it as barebone as possible. The problem boils down to τend not being defined due to the error in Event Handling. I plot the functions at the end with the surface where the integration should stop in black. Note that in this simpler example, the interpolation functions for the coordinates are real but it still doesn't solve the problem.

coords = {t,r,θ,ϕ};
n=Length[coords];
rplus = M+Sqrt[M^2-a^2];
tt=2 M r/ρ-1;
rr=ρ/Δ;
θθ=ρ;
ϕϕ=(Δ+(2 M r (r^2+a^2))/ρ)Sin[θ]^2;
tϕ=-4 a M r Sin[θ]^2/ρ;
metric={{tt,0,0,tϕ},{0,rr,0,0},{0,0,θθ,0},{tϕ,0,0,ϕϕ}};
metric//MatrixForm;
inversemetric=Simplify[Inverse[metric]];
inversemetric //MatrixForm;

a=0;M=.;
ρ=r^2+a^2 Cos[θ]^2;
Δ=r^2 - 2 M r+ a^2;
christoffel:=christoffel=Simplify[Table[(1/2)*Sum[(inversemetric[[i,s]])*
(D[metric[[s,j]],coords[[k]] ]+
D[metric[[s,k]],coords[[j]] ]-D[metric[[j,k]],coords[[s]] ]),{s,1,n}],   {i,1,n},{j,1,n},{k,1,n}] ]
listchristoffel:=Table[If[UnsameQ[christoffel[[i,j,k]],0],   {ToString[Γ[i,j,k]],christoffel[[i,j,k]]}] ,{i,1,n},{j,1,n},{k,1,j}]
geodesic:=geodesic=Simplify[Table[-Sum[christoffel[[i,j,k]]coords[[j]]'      coords[[k]]',{j,1,n},
{k,1,n}],{i,1,n}]]
listgeodesic:=Table[{"d/dτ"ToString[coords[[i]]'],"=",geodesic[[i]]},i,1,n}]
geodesic:=geodesic=Simplify[Table[-Sum[christoffel[[i,j,k]]coords[[j]]' coords[[k]]',{j,1,n},
{k,1,n}],{i,1,n}]]
TableForm[Partition[DeleteCases[Flatten[listchristoffel],Null],2],TableSpacing->{2,2}]
TableForm[Partition[DeleteCases[Flatten[listgeodesic],Null],2],TableSpacing->{2,2}]

uinvar=0;M=1; τval=20; rplus = M + Sqrt[M^2-a^2];
ivs = {-0.15,0.0,0.001}; ics = {0,5,π/3,π/4};
ivs=Join[{χ},ivs];
tmp=metric;
tmp=tmp/.Table[coords[[i]]->ics[[i]],{i,0,n}];
tmp=ivs.(tmp.ivs);
χslv=Solve[tmp==uinvar,χ];
ivs[[1]]=Last[χ/.χslv];
deq=Table[coords[[i]]''[τ]
          ==(geodesic[[i]]/.Join[
                                 Table[coords[[i]]'->coords[[i]]'[τ],{i,1,n}],
                                 Table[coords[[i]]->coords[[i]][τ],{i,1,n}]]),{i,1,n}];
deq=Join[deq,Table[coords[[i]]'[0]==ivs[[i]],{i,1,n}],
             Table[coords[[i]][0]==ics[[i]],{i,1,n}]];
soln=NDSolve[deq,coords,{τ,0,τval},
             Method->{"EventLocator",
                      "Event":>Re[coords[[2]][τ]]=2,
                      "EventAction":>Throw[τend=τ,"StopIntegration"]}];

Edit 2: I thought I had solved this by adding a

[[1]]

behind the

coords[[2]][τ]

in the "Event" but this wasn't correct. It does get rid of the error, but

coords[[2]][τ][[1]]

seems to call the argument of the function... As can be seen by calling

coords[[2]][5][[1]]

for example.

I think I have narrowed the problem down to being about how I'm calling the coords list, but I seem to have made such a mess of it that I can't figure it out.

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Dr. belisarius Mar 20 '15 at 16:56
  • $\begingroup$ Try for example "Event" :> ((Re[r[\[Tau]]] - 2) == 0) .... $\endgroup$ – Dr. belisarius Mar 20 '15 at 17:02
  • $\begingroup$ Now the error message disappears but the == won't suffice as I doubt that difference will ever be exactly zero. However, changing it into <= should do the trick. Thanks! It's however not recognizing it passing that particular radius... In fact, this Event doesn't catch anything. Even if I put the 2 to 100. $\endgroup$ – helgi Mar 20 '15 at 17:16
  • $\begingroup$ After a bit of playing around with this, the trick seems to have been with using r[[Tau]] rather than coords itself. Thank you so much! For future readers, the Event line now reads Method -> {"EventLocator", "Event" :> (r[[Tau]] <= 2), "EventAction" :> Throw[[Tau]end = [Tau], "StopIntegration"]}]; $\endgroup$ – helgi Mar 20 '15 at 17:23
  • $\begingroup$ Please answer your own question for posterity. Thanks! $\endgroup$ – Dr. belisarius Mar 20 '15 at 17:29
1
$\begingroup$

The solution was found in referencing the actual coordinate

r[\[Tau]]

rather than the its placement in the list coords.

$\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.