Twee pieken of toch maar één?

Moderator: physicalattraction

Reageer
Gebruikersavatar
Berichten: 6.716

Re: Twee pieken of toch maar één?

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.

Gebruikersavatar
Berichten: 6.716

Re: Twee pieken of toch maar één?

Hier de bijbehorende plots:

screenshot.png

Berichten: 28

Re: Twee pieken of toch maar één?

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

Technicus
Berichten: 904

Re: Twee pieken of toch maar één?

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.

Gebruikersavatar
Berichten: 6.716

Re: Twee pieken of toch maar één?

wnvl1 schreef: do 22 jul 2021, 18:20 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.

Gebruikersavatar
Berichten: 6.716

Re: Twee pieken of toch maar één?

CoenCo schreef: do 22 jul 2021, 18:38 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:

tien-maal.png

Gebruikersavatar
Berichten: 6.716

Re: Twee pieken of toch maar één?

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:
vergeleken.png

Gebruikersavatar
Berichten: 6.716

Re: Twee pieken of toch maar één?

Even terug in de tijd:
Professor Puntje schreef: vr 11 jun 2021, 23:39 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 = Rzon en r = √(x2 + Rzon2). 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...

Berichten: 28

Re: Twee pieken of toch maar één?

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.

Berichten: 28

Re: Twee pieken of toch maar één?

Afbeelding

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

Berichten: 28

Re: Twee pieken of toch maar één?

Afbeelding

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.

Berichten: 28

Re: Twee pieken of toch maar één?

wnvl1 schreef: vr 23 jul 2021, 02:55 Afbeelding

Dus dit is het coördinatenstelsel. Maar is het zeker dat Solomon Antoniou dat ook gebruikt?
Nee mijn y en x moeten omgekeerd staan.

Gebruikersavatar
Berichten: 6.716

Re: Twee pieken of toch maar één?

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.

Gebruikersavatar
Berichten: 6.716

Re: Twee pieken of toch maar één?

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

Gebruikersavatar
Berichten: 6.716

Re: Twee pieken of toch maar één?

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?

Reageer