Twee pieken of toch maar één?

[Professor Puntje]

Onder voorbehoud (hier zitten wellicht nog fouten in):

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ellipj
from scipy.special import ellipk

fig, ax= plt.subplots(4,1,figsize=(4, 6))



# Berekening lichtbaan

c = 3e8
r0 = 7e8
rs = 2.95e3

e1 = (r0 - rs + np.sqrt((r0 - rs)*(r0 + 3*rs)))/(2*rs*r0)
e2 = 1/(r0)
e3 = (r0 - rs - np.sqrt((r0 - rs)*(r0 + 3*rs)))/(2*rs*r0)
tau = np.sqrt((rs*(e1 - e3))/4)
h = np.sqrt((e2 - e3)/(e1 - e3))
m = h*h
sigma = -tau*(np.pi/2) + ellipk(m)

t = np.linspace(np.pi/2-1.4,np.pi/2+1.4,1000)

sn,_,_,_=ellipj(tau*t + sigma , m)

x = (1/(e3 + (e2 - e3)*sn**2))*np.cos(t)
y = (1/(e3 + (e2 - e3)*sn**2))*np.sin(t)



# Berekening diverse afgeleiden

dx=np.diff(x)
dy=np.diff(y)
dydx=dy/dx

r=np.sqrt(x**2 + y**2)
S=(x[:-1] + y[:-1]*dydx)/r[:-1]
T=(S*(x[:-1]/r[:-1]) - 1)/y[:-1]
dx2dt2=(( (1 - (rs/r[:-1]))**2 )/( S**2 + (r[:-1]**2)*(1 - (rs/r[:-1]))*(T**2) ))*(c**2)

ddydx=np.diff(dydx)
ddydxdx=ddydx/dx[:-1]



# Plots

ax[0].plot(x,y,c=(1,0.2,0.5),lw=1)
ax[0].title.set_text('y = f(x)')

ax[1].plot(x[:-1],dydx,c=(1,0.2,0.5),lw=1)
ax[1].title.set_text('dy/dx')

ax[2].plot(x[:-1],dx2dt2,c=(1,0.2,0.5),lw=1)
ax[2].title.set_text('dx2/dt2')

ax[3].plot(x[:-2],ddydxdx,c=(1,0.2,0.5),lw=1)
ax[3].title.set_text('d2y/dx2')

plt.tight_layout()
plt.show()

Graag commentaar van ervaren Python gebruikers.

[Professor Puntje]

Hier de bijbehorende plots:

[wnvl1]

Is er een reden dat je gebruik maakt in (9.19) en niet (9.11) uit dat artikel?

[CoenCo]

Ik heb niet alle formule's gechecked, maar zag wel deze:
Een van de dingen waar nog wat onnauwkeurigheid in kan zitten is het gebruik van np.diff()
Dit is een backward difference, en daarmee niet de meest nauwkeurige methode. Eigenlijk bepaal je steeds de afgeleide ongeveer midden tussen 2 punten in, i.p.v. op het punt zelf.
scipy heeft ook een central difference, die mogelijk nauwkeuriger is.

Voor uitleg over het verschil, zie bijvoorbeeld:
https://pythonnumericalmethods.berkeley ... tives.html

En voor de functies:
https://numpy.org/doc/stable/reference/ ... .diff.html
https://docs.scipy.org/doc/scipy/refere ... ative.html

Overigens is vrij makkelijk te testen of dit invloed heeft, zonder een rewrite van je code: Door het verdubbelen van het aantal stappen in je linspace (nu 1000 voor t), halveer je de fout in de numerieke differentiatie.

[Professor Puntje]

Is er een reden dat je gebruik maakt in (9.19) en niet (9.11) uit dat artikel?

Ja - voor de Weierstrass P-functie heb ik geen geschikte Linux programmatuur gevonden.

[Professor Puntje]

Overigens is vrij makkelijk te testen of dit invloed heeft, zonder een rewrite van je code: Door het verdubbelen van het aantal stappen in je linspace (nu 1000 voor t), halveer je de fout in de numerieke differentiatie.

Met een extra nulletje (dus 10000 voor t) krijgen we:

[Professor Puntje]

Ik heb nog even een plot van formule 6 aan het progje toegevoegd:

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ellipj
from scipy.special import ellipk

fig, ax= plt.subplots(5,1,figsize=(4, 6))



# Berekening lichtbaan

c = 3e8
r0 = 7e8
rs = 2.95e3

e1 = (r0 - rs + np.sqrt((r0 - rs)*(r0 + 3*rs)))/(2*rs*r0)
e2 = 1/(r0)
e3 = (r0 - rs - np.sqrt((r0 - rs)*(r0 + 3*rs)))/(2*rs*r0)
tau = np.sqrt((rs*(e1 - e3))/4)
h = np.sqrt((e2 - e3)/(e1 - e3))
m = h*h
sigma = -tau*(np.pi/2) + ellipk(m)

t = np.linspace(np.pi/2-1.4,np.pi/2+1.4,1000)

sn,_,_,_=ellipj(tau*t + sigma , m)

x = (1/(e3 + (e2 - e3)*sn**2))*np.cos(t)
y = (1/(e3 + (e2 - e3)*sn**2))*np.sin(t)



# Berekening diverse afgeleiden

dx=np.diff(x)
dy=np.diff(y)
dydx=dy/dx

r=np.sqrt(x**2 + y**2)
S=(x[:-1] + y[:-1]*dydx)/r[:-1]
T=(S*(x[:-1]/r[:-1]) - 1)/y[:-1]
dx2dt2=(( (1 - (rs/r[:-1]))**2 )/( S**2 + (r[:-1]**2)*(1 - (rs/r[:-1]))*(T**2) ))*(c**2)

form6=(( ( 1 - (rs/r) )**2 )/( (rs/r)*( (x**2)/(r**2) ) + 1 - rs/r ) )*(c**2)

ddydx=np.diff(dydx)
ddydxdx=ddydx/dx[:-1]



# Plots

ax[0].plot(x,y,c=(1,0.2,0.5),lw=1)
ax[0].title.set_text('y = f(x)')

ax[1].plot(x[:-1],dydx,c=(1,0.2,0.5),lw=1)
ax[1].title.set_text('dy/dx')

ax[2].plot(x[:-1],dx2dt2,c=(1,0.2,0.5),lw=1)
ax[2].title.set_text('dx2/dt2 (exact)')

ax[3].plot(x,form6,c=(1,0.2,0.5),lw=1)
ax[3].title.set_text('dx2/dt2 (form. 6)')

ax[4].plot(x[:-2],ddydxdx,c=(1,0.2,0.5),lw=1)
ax[4].title.set_text('d2y/dx2')

plt.tight_layout()
plt.show()

Dat geeft:

[Professor Puntje]

Even terug in de tijd:

Ik wil eens zien of we wat van de toegepaste benaderingen in de afleiding kunnen omzeilen. Terug naar formules (6) en (7):

\(\)

\( \frac{\mathrm{d}x^2}{\mathrm{d}t^2} = \frac{(1 - \frac{r_s}{r})^2}{ \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} } \cdot c^2 \,\,\,\,\,\,\,\, (6)\)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{\frac{\partial}{\partial R} \left ( \frac{\mathrm{d} x}{\mathrm{d}t} \right ) }{\frac{\mathrm{d}x}{\mathrm{d}t}} \,\,\,\,\,\,\,\, (7) \)

\(\)

Dat geeft:

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{\frac{\partial}{\partial R} \left ( \frac{\mathrm{d} x}{c \, \mathrm{d}t} \right ) }{\frac{\mathrm{d}x}{c \, \mathrm{d}t}} \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{\partial}{\partial R} \ln(\frac{\mathrm{d}x}{c \, \mathrm{d}t}) \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{\partial}{\partial R} \ln \left (\sqrt{\frac{\mathrm{d}x^2}{c^2 \mathrm{d}t^2}} \right ) \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{1}{2} \frac{\partial}{\partial R} \ln(\frac{\mathrm{d}x^2}{c^2 \mathrm{d}t^2}) \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{1}{2} \frac{\partial}{\partial R} \ln \left (\frac{(1 - \frac{r_s}{r})^2}{ \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} } \right ) \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{1}{2} \frac{\partial}{\partial R} \left \{ \ln \left ((1 - \frac{r_s}{r})^2 \right ) \, - \, \ln \left ( \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} \right )\right \} \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{1}{2} \frac{\partial}{\partial R} \ln \left ((1 - \frac{r_s}{r})^2 \right ) \, - \, \frac{1}{2} \frac{\partial}{\partial R} \ln \left ( \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} \right ) \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{\partial}{\partial R} \ln (1 - \frac{r_s}{r}) \, - \, \frac{1}{2} \frac{\partial}{\partial R} \ln \left ( \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} \right ) \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{\frac{\partial}{\partial R} (1 - \frac{r_s}{r})}{ 1 - \frac{r_s}{r} } \, - \, \frac{1}{2} \frac{ \frac{\partial}{\partial R} \left ( \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} \right ) }{ \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} } \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{- r_s \frac{\partial}{\partial R} \frac{1}{r}}{ 1 - \frac{r_s}{r} } \, - \, \frac{1}{2} \frac{r_s x^2 \frac{ \partial}{\partial R} \frac{1}{r^3} - r_s \frac{\partial}{\partial R} \frac{1}{r} }{ \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} } \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \frac{ r_s \frac{R}{r^3}}{ 1 - \frac{r_s}{r} } \, - \, \frac{1}{2} \frac{- 3 r_s x^2 \frac{R}{r^5} + r_s \frac{R}{r^3} }{ \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} } \)

\(\)

\( \frac{\mathrm{d}\varphi}{\mathrm{d}x} = \left ( \frac{1}{ 1 - \frac{r_s}{r} } \, - \, \frac{1}{2} \frac{- 3 \frac{x^2}{r^2} + 1 }{ \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} } \right ) \cdot r_s \frac{R}{r^3} \,\,\,\,\,\,\, (16) \)

\(\)

De vraag is nu wat voor curve dit oplevert voor R = R_zon en r = √(x² + R_zon²). Zitten ook hier die twee pieken al in...?.

En toen werden in een plot van (16) ook twee pieken gevonden. De route van formule (6) naar formule (16) loopt afgezien van de toepassing van Huygens' principe echter zonder verdere benaderingen. Als de grafiek van dφ/dx als functie van x inderdaad maar één piek mag bevatten zou Huygens' principe de twee pieken van \( \frac{dx^2}{dt^2} \) als functie van x ongedaan moeten maken en in één piek moeten omzetten, iets wat kennelijk niet gebeurt...

[wnvl1]

In de paper van Mathpages gaat het over licht dat van oneindig ver komt en wordt afgebogen door de zon als ik het juist begrijp. In de paper van Solomon Antoniou gaat het over licht dat in een orbitaal rond de zon draait??? Zijn beiden dan wel zomaar te vergelijken? Het is gewoon een bedenking, zeg niet dat het niet zo is.

Waaruit kan je opmaken dat
X=r*cos(t)
Y=r*sint(t)
en niet omgekeerd? Het is spijtig dat in de paper van Solomon Antoniou geen enkele figuur met coördinaten staat.

[wnvl1]

Dus dit is het coördinatenstelsel. Maar is het zeker dat Solomon Antoniou dat ook gebruikt?

[wnvl1]

Het licht gaat niet in een orbitaal rond de zon, maar toch plaats de auteur het in een sectie met als titel orbitalen... Dit zou wel kunnen bij een zwart gat. Al lijkt het mij voor de berekening waarschijnlijk niet veel uit te maken.

[wnvl1]

Dus dit is het coördinatenstelsel. Maar is het zeker dat Solomon Antoniou dat ook gebruikt?

Nee mijn y en x moeten omgekeerd staan.

[Professor Puntje]

De auteur van het artikel behandelt zowel gesloten lichtbanen als niet-gesloten lichtbanen. De formule die ik gebruik beschrijft de niet-gesloten lichtbanen in de Schwarzschild coördinaten φ en r. De pseudo-cartesische coördinaten x en y bereken ik hier zelf uit φ en r met behulp van de transformatie:

\(\)

\( x = r \cos(\alpha) \\ y = r \sin(\alpha) \)

\(\)

Die x en y komen niet uit het artikel.

[Professor Puntje]

Bedenk ook dat x en y de werkelijke afstanden vertekend weergeven, we hebben hier immers met gekromde ruimtetijd van doen.

[Professor Puntje]

Dit topic heeft al zoveel verrassende wendingen gekend dat ik wat huiverig ben om nu weer een voorlopige conclusie te trekken, maar zonder dat komen we ook niet verder. Dus vooruit dan maar:

Het ziet ernaar uit dat de aanpak van MathPages (in de hier gereconstrueerde vorm) ondanks de vele benaderingen een juist resultaat voor \( \frac{dx^2}{dt^2} \) als functie van x oplevert en dat er in de grafiek van die functie twee (neerwaartse) pieken voorkomen.

Graag commentaar:
- Mogen we dat nu inderdaad als voorlopige conclusie beschouwen?
- Of moet er daarvoor eerst nog iets extra's onderzocht worden?