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.
"Event" :> ((Re[r[\[Tau]]] - 2) == 0).... $\endgroup$