Part 1
I am trying to integrate Hamiltons equations for the Schwarzschild geometry using NDSolve, these equations must all be integrated simultaneously and are nice ODE's.
A similar task is done in this paper but for a wormhole geometry, Interstellar Wormholes look at appendix page 12 for the equations.
Via equivalent calculation you can get 5 Hamiltons equations for the Schwarzschild geometry calculation. I want to be able to integrate my analogous equations like is spoken about in the paper for their equations (at the end of page 12 and beginning of page 13).
For some background:
The equivalent Super-Hamiltonian for my case is
$\qquad H=\frac{1}{2}[\frac{1}{1-r_s/r}p_t^2-(1-r_s/r)p_r^2-\frac{1}{r^2}p_\theta^2-\frac{1}{r^2sin^2(\theta)}p_\phi^2]$
where $r_s=2M$ and $M$ is some arbitrary mass.
- Thus the resulting "Hamiltons equations" are given by the following where $H$ is defined above. Note the $p_t$ and $p_\phi$ equations are 0 since the Hamiltonian is independent of those variables. This results in trival solutions for them; i.e., $p_t=E$ and $p_\phi=b$, where $E$ and $b$ are constants:
$\qquad dt/d\lambda=\frac{\partial H}{\partial p_t};$
$\qquad dr/d\lambda=\frac{\partial H}{\partial p_r};$
$\qquad d\theta/d\lambda=\frac{\partial H}{\partial p_\theta};$
$\qquad d\phi/d\lambda=\frac{\partial H}{\partial p_\phi};$
$\qquad dp_t/d\lambda=-\frac{\partial H}{\partial t};$
$\qquad dp_r/d\lambda=-\frac{\partial H}{\partial r};$
$\qquad dp_\theta/d\lambda=-\frac{\partial H}{\partial \theta};$
$\qquad dp_\phi/d\lambda=-\frac{\partial H}{\partial \phi};$
Once these are calculated it emerges that some are easily integrated due to symmetries and conservation laws.
I am not quite sure how to input these equations correctly along with the necessary initial conditions and constants into NDSolve (I am typing the equations manually into Mathematica but not getting very far, not sure if its how I am typing the equations in which is the problem).
This is what I have tried doing so far.
Defining a list of the differential equations (from Hamiltons equations above) to be put into Mathematica so they can be integrated simultaneously (ham as in Hamiltonian).
ham = {dt/dλ + E/(1 - (2 M)/r) == 0, dr/dλ - (1 - (2 M)/r) Subscript[p, r] == 0, dθ/dλ - Subscript[P, θ]/r^2 == 0, dϕ/dλ - b/(r^2 (Sin^2)[θ]) == 0, Subscript[dp, r]/dλ + M/r^2 (E^2/(1 - (2 M)/r)^2 + (Subscript[p, r])^2) - B^2/r^3 == 0, Subscript[dp, θ]/dλ - (b^2 Cos[θ])/(r^2 (Sin^3)[θ]) == 0}Defining a list of the initial conditions for the integration of these equations since they must be solved numerically, since they are 1st order DE they only need 1 piece of info for each one ie initial function values. This might be wrong, just choose random initial conditions.
haminital = {t[0] = 3, r[0] = 10, θ[0] = 2, ϕ[0] = Pi/2, Subscript[p, r][0] = 4, Subscript[p, θ][0] = 10}Putting it all into NDSolve along with constants E, b and B seperatley to try and solve for our functions.
NDSolve[{ham, haminital}, {t, r, θ, ϕ, Subscript[p, r], Subscript[p, θ]}, ?]
Not sure what to put in question mark spot.
All in all I am not even sure if I have got the correct initial conditions or if I am integrating the right equations or not.
Part 2
I am now attempting to do this using the correct initial conditions set up as shown on page 13 (Procedure for generating map) of Interstellar Wormholes. This is so I can generate maps for trajectories.
Can someone help piece this all together?
The first thing it says (step 1) is choose a camera location (I am keeping it as variables so I can input specific values for them in later.
camloc = {rc, θc, ϕc};The second thing it says (step 2) is set up cartesian coordinates entered on the cameras location we defined in "step 1". Not exactly sure how to create cartesian coordinates that are centred on my
{rc, θc, ϕc}.The third thing it says (step 3) is to setup spherical polar coordinates from those cartesian coordinates which are centred on our choice of camera location "step 2". Once these are set up choose some point in the cameras sky
{θcs, ϕcs}This isn't too bad I guess once you have set up your standard cartesian coordinate since you can use coordinate transformation equations to get them interms of spherical polar, my guess is something like this:
Nx = Sin[θcs] Cos[ϕcs]; Ny = Sin[θcs] Sin[ϕcs]; Nz = Cos[θcs]; camsky = {Nx, Ny, Nz}Note: Nx,Ny,Nz are unit vector components so we get r=1 and due to spherical symmetry. Also $\theta cs$ and $\phi cs$ just represent the angles seen from the cameras perspective (Coordinates which should have been setup in step 2) and in-order to not confuse with the $\theta$ and $\phi$ of the metric coordinate which we will us later.
These are then related to the components of a unit vector for a light ray (nr, n$\phi$, n$\theta$) by:
nr = -Nx; nϕ = -Ny; nθ = Nz;Then it says is compute the canonical momenta from the following set of relations (These have been altered from the text since we are using a different metric).
pr = (1 - 2 M/r)^-1 nr; pθ = (1 - 2 M/r)^-1/2 rnθ; pϕ = (1 - 2 M/r)^-1/2 rSin[θ] nϕ;Note: These are the same canonical momenta that I and Alex defined when inputting the initial conditions into
NDSolve. Also, theθused in the final condition is not the same as theθ cswe defined earlier it is instead the coordinate $\theta$ I spoke about.Also, since the actual constants $b$ and $B^2$ are related to the canonical momenta "we just defined above", we get:
b = pϕ; B^2 = pθ^2 + pϕ^2/Sin^2[θ];Note: Now all of the canonical momenta (pr.. etc) are found and since the constants b and B^2 depend on the canonical momenta those are also defined. The other constants E0 and M can be chosen semi-arbitrarily so don't matter too much.
This is the final thing (step 6) where we now take everything we have just calculated and bring it together and use it as initial conditions.
It says take as initial conditions:
1) For the positions: since we want at $\lambda=0$ the light ray starts at the cameras position (defined in step 1)
{r[0]==rc, θ[0]==θc , ϕ[0]==ϕc}2) For canonical momenta $(Pr,P\theta,P\phi)$: since we want at $\lambda=0$ the momenta that we defined in step 5
{pr[0] == (1 - 2 M/r)^-1 nr, pθ[0] == (1 - 2 M/r)^-1/2 rnθ, pϕ[0] == (1 - 2 M/r)^-1/2 rSin[θ] nϕ}3) For constants $(b,B^2)$:
{b == pϕ[0], B^2 == pθ[0]^2 + pϕ[0]^2/Sin^2[θ]Note: I think I may need to put all of the $r,\theta,\phi,Pr$.. etc as functions of lambda for it to work properly (attempted below). This is then all integrated in the range $\lambda$ = some large negative number to $\lambda=0$ using
NDSolve
This is what I potentially think it could look like (However I didn't quite get step 2 so it is wrong but hopefully correctable).
camloc = {rc, θc, ϕc};
Nx = Sin[θcs]*Cos[ϕcs];
Ny = Sin[θcs]*Sin[ϕcs];
Nz = Cos[θcs];
camsky = {Nx, Ny, Nz}
nr = -Nx;
nϕ = -Ny;
nθ = Nz;
pr[λ_] = (1 - 2 M/r[λ])^-1 nr;
Pθ[λ_] = (1 - 2 M/r[λ])^-1/2 r[λ] nθ;
pϕ[λ_] = (1 - 2 M/r[λ])^-1/2 r[λ] Sin[θ[λ]] nϕ;
b = pϕ[λ];
B^2 = pθ[λ]^2 + pϕ[λ]^2/Sin^2[θ[λ]];
M = 1; E0 = .1;
ham =
{t'[λ] + E0/(1 - (2 M)/r[λ]) == 0,
r'[λ] - (1 - (2 M)/r[λ]) pr[λ] == 0,
θ'[λ] - Pθ[λ]/r[λ]^2 == 0,
ϕ'[λ] - b/(r[λ]*Sin[θ[λ]])^2 == 0,
pr'[λ] + M/r[λ]^2 (E0^2/(1 - (2 M)/r[λ])^2 + pr[λ]^2) - B^2/r[λ]^3 == 0,
Pθ'[λ] - (b^2*Cos[θ[λ]])/(r[λ]^2 Sin[θ[λ]]^3) == 0};
haminital =
{t[0] == 3, r[0] == rc, θ[0] == θc, ϕ[0] == ϕc, pr[0] == 4, Pθ[0] == 10};
sol =
NDSolve[{ham, haminital}, {t, r, θ, ϕ, pr, Pθ}, {λ, -100000, 0}]
I need help debugging it (filling in steps 2 and general). Feel free to put in values for $rc,\theta c,\phi c$ along with $rcs,\theta cs,\phi cs$ to get some plots.
Part 3
I am trying to now plot these numerical solutions plotted using ParametricPlot3D where $\lambda$ is the parameter.
Not quite sure why this doesn't work
gx[λ_] == (r[λ_])*Sin[θ[λ_]]*
Cos[ϕ[λ_]];
gy[λ_] == (r[λ_])*Sin[θ[λ_]]*
Sin[ϕ[λ_]];
gz[λ_] == (r[λ_])*Cos[θ[λ_]];
ParametricPlot3D[{gx[λ_], gy[λ_],
gz[λ_]}, {λ, -100, 0}]
Where the $r, \theta$ and $\phi$ functions are the solutions we got from the numerical integration. Also, we have converted our $r, \theta$ and $\phi$ to $x,y,x$ by spherical coordinate transformations (since Mathematica only recognises cartesian values).
Part 4 (Currently on)
On page 13 of the paper bullet point 3 "Implementing map" It says and used that map, numerical table of final angle values as functions of our initial camera angles in sky $(\phi cs, \phi[\lambda])$".
Having already generated a single photon ray (Our solutions to the integration of angles etc for 1 set of initial conditions) how do I create a table of these solutions for a range of angles all in 1 go?
This is my attempt to list out the values of $\phi cs$ using a four-loop along with the our value for $\phi[\lambda]$ for the end of our integration range (-100).
numericalmap = {};
numericalmap;
For[i = 1, i < 50, i++,
\[Phi]csgen = {1 + 0.1 i} \[Phi]cs;
\[Phi]2[\[Lambda]_] = \[Phi][\[Lambda]] /. sol[[1]];
numericalmap = Append[numericalmap, {\[Phi]csgen, \[Phi]2[-100]}];
]
numericalmap
The problem is that I can't get my numerical value of $\phi [-100]$ to change as the four-loop puts in different $\phi cs$ values.
My whole thing looks like, could someone debug please
camloc = {rc, θc, ϕc};
Nx = Sin[θcs]*Cos[ϕcsgen];
Ny = Sin[θcs]*Sin[ϕcsgen];
Nz = Cos[θcs];
camsky = {Nx, Ny, Nz};
nr = -Nx;
nϕ = -Ny;
nθ = Nz;
(*pr[λ_]=(1-2 M/r[λ])^-1 nr;
Pθ[λ_]=(1-2 M/r[λ])^-1/2 r[λ]*nθ;
pϕ[λ_]=(1-2 M/r[λ])^-1/2 r[λ]*Sin[\
θ[λ]] nϕ;*)
(*b=pϕ[λ];
B^2=pθ[λ]^2+pϕ[λ]^2/Sin^2[θ[\
λ]];*)
b = rcs*Sin[θcs]*nϕ;
B2 = rcs^2*(nϕ^2 + nθ^2);
M = 1; E0 = .1; {rcs, θcs, ϕcs} = {1000, Pi/2, Pi};
ham = {t'[λ] + E0/(1 - (2 M)/r[λ]) == 0,
r'[λ] - (1 - (2 M)/r[λ]) pr[λ] ==
0, θ'[λ] - Pθ[λ]/r[λ]^2 ==
0, ϕ'[λ] -
b/(r[λ]*Sin[θ[λ]])^2 == 0,
pr'[λ] +
M/r[λ]^2 (E0^2/(1 - (2 M)/r[λ])^2 +
pr[λ]^2) - B2/r[λ]^3 == 0,
Pθ'[λ] - (b^2*
Cos[θ[λ]])/(r[λ]^2*
Sin[θ[λ]]^3) == 0};
haminital = {t[0] == 3,
r[0] == rcs, θ[0] == θcs, ϕ[0] == ϕcs,
pr[0] == 4, Pθ[0] == 10};
sol = NDSolve[{ham, haminital}, {t, r, θ, ϕ, pr,
Pθ}, {λ, 0, -100}]
numericalmap = {};
numericalmap;
For[i = 1, i < 50, i++,
\[Phi]csgen = {1 + 0.1 i} \[Phi]cs;
\[Phi]2[\[Lambda]_] = \[Phi][\[Lambda]] /. sol;
numericalmap = Append[numericalmap, {\[Phi]csgen, \[Phi]2[-100]}];
]
numericalmap
