0
$\begingroup$

These are the equations for Poincaré sections:

    M = 0.5; a = 0.2; s = 0.1;M = 1/2; Δ = r^2 - 2*M*r[t] + a^2;
    W = r[t]^2 + a^2*Cos[θ[t]]^2;
    h = s*M^3*r[t]/W^2;
    gtt = -Simplify[(1 - 2*M*r[t]/W)*(1 + h)];
    grr = Simplify[W (1 + h)/(Δ + a^2*Sin[θ[t]]^2*h)];
    gss = W;
    gff = Simplify[(r[t]^2 + a^2 + 2 a^2*M*r[t]*Sin[θ[t]]^2/W)*
    Sin[θ[t]]^2 + h*a^2*Sin[θ[t]]^4*(1 + 2*M*r[t]/W)]
    gtf = -2*a*M*r[t]*Sin[θ[t]]^2*(1 + h)/W
    L = 1/2 (gtt*en^2 + grr*D[r[t], t]^2 + gss*D[θ[t], t]^2 + 
    gff*ln^2 + 2 gtf*en*ln);

I would like to plot the Poincaré section for collection of (r[t], Derivative[1][r][t]), I wrote:

   date = Block[{en = 0.95, ln = 1.5}, 
   Reap[NDSolve[{grr*(r^′′)[t] == 
   D[L, r[t]] - D[grr, r[t]]*Derivative[1][r][t]^2 - 
    D[grr, θ[t]]*Derivative[1][r][t]*
     Derivative[1][θ][t], 
  gss*(θ^′′)[t] == 
   D[L, θ[t]] - 
    D[gss, θ[t]]*Derivative[1][θ][t]^2 - 
    D[gss, r[t]]*Derivative[1][r][t]*Derivative[1][θ][t], 
  r[0] == 8, θ[0] == Pi/2,
  WhenEvent[Derivative[1][r][0] == 0, 
   Sow[{r[t], Derivative[1][r][t]}]]}, {}, {t, 100000}, 
 MaxSteps -> ∞]]][[-1, 1]]; ListPlot[data, 
ImageSize -> Medium, PlotRange -> {{0, 20}, All}, 
 PlotStyle -> PointSize[0.0025]]

But it does not work.

$\endgroup$
2
$\begingroup$

The code in the question contain (r^′′)[t] and (θ^′′)[t], which are, I presume, transcription errors. They should be r''[t] and θ''[t]. With these corrections, the code contains two errors. First, Δ should be given by

Δ = r[t]^2 - 2*M*r[t] + a^2

Second, values for r'[0]and θ'[0] must be provided to NDSolve. For all values of r'[0]and θ'[0] that I have tried, r[t] grows exponentially, and at most one value of {r[t], r'[t]} is returned by WhenEvent. So, for instance, with r'[0] == -10, θ'[0] == -2, one obtains (with {} replaced by {r[t], θ[t]} in NDSolve and [[-1, 1]] removed),

LogPlot[Evaluate[{r[t], θ[t]} /. date[[1, 1]]], {t, 0, 10}]

enter image description here

shows that r[t] grows exponentially and θ[t] becomes constant after only a single oscillation. Thus,

date[[-1, 1]]
(* {{6.96449, 1.87522*10^-12}} *)

Perhaps a better choice for r'[0]and θ'[0] would yield an oscillatory solution, but finding it is not trivial.

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