Formule (9) en Python

[Professor Puntje]

Beste Professor,

Ik ben al de hele ochtend uren jouw methode aan het bestuderen. Het werkt echt makkelijker als je zelf ook leert een plotje te maken en jouw resultaten te controleren. Nu moet ik mijn onkunde toepassen op jouw afleidingen om daadwerkelijk bevestiging te krijgen of jouw formule klopt.

Mijn analyse jouw aanpak:
$$\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} \,\,\,\,\,\,\, (9)$$

Mijn input:

Code: Selecteer alles

 phi=(1/(1-Rs/r)-0.5*(-3*(x**2/r**2)+1)/((x**2/r**2) *Rs/r  +1 -Rs/r))*Rs*R/r**3

Code: Selecteer alles

#Open pyplot in separate interactive window
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt5')

#https://www.mathpages.com/rr/s8-09/8-09.htm
#https://www.mathpages.com/rr/s6-03/6-03.htm

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams

rcParams['axes.titlepad'] = 20 

widths = [10,10]
heights = [10]

fig= plt.figure(figsize=(20,10))

gs=fig.add_gridspec(1,2,width_ratios=widths, height_ratios=heights)

ax1=fig.add_subplot(gs[0,0])
ax2=fig.add_subplot(gs[0,1])


#Physical constants
M=1.989e30
G=6.67408e-11
c=3e8
R=696340000

#schwarzschild Radius
Rs=4*M*G/c**2

#Function c(r)
def fdphidr1(r,x):
    
    phi=(1/(1-Rs/r)-0.5*(-3*(x**2/r**2)+1)/((x**2/r**2) *Rs/r  +1 -Rs/r))*Rs*R/r**3
  
    return phi


y=R
x=np.linspace(-10,10,10000)
x=x*R
r=np.sqrt(y**2+x**2)

#Angular Distribution
dphidr=fdphidr1(r,x)

ax1.plot(x/R,dphidr,color="red", linewidth=0.5, label=r"(9)")
ax1.set_xlabel(r'$x$',fontsize=15)
ax1.set_ylabel(r'$d \phi /dx$',fontsize=15)

#Integrated deflection agngle
deflection=np.cumsum(dphidr)*(x[10]-x[9])
deflection=np.degrees(deflection)*3600

ax2.plot(x/R,deflection,color="red", linewidth=0.5, label=r"(9)")
ax2.set_xlabel(r'$x$',fontsize=15)
ax2.set_ylabel(r'$\int \phi dx$',fontsize=15)


ax1.legend(loc="upper right")
ax2.legend(loc="upper right")

Puntje.png

Mijn obervatie:

• Er zijn twee pieken.

• Totale deflectiehoek krijg ik niet kloppende.

Ik ben jouw afleiding nagelopen.

Volgens mij zou binnen diff. deze formule van jouw \(c(r)\) dienen te zijn. Ik snap niet waar die \(\ln\) vandaan komt (deze formule krijg ik kloppend als ik \(\ln\) door \(\sqrt{}\) vervang).
$$\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 )$$
Ergens diff. je \(1/R\) naar \(\ln\) mag dat? \(dx\) is niet liniear met \(dR\)?
Ik snap niet wat ik fout doe en snap het niet. Totaal overmeesterd door mijn onkunde.

Volgens mijn simpele ziel is en mijn kennis van GR is \(dt^2\) \((dt)^2\):
$$\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)$$
Is dit niet:
$$\frac{\mathrm{d}x}{\mathrm{d}t} = \sqrt{\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)$$
Mijn simpele ziel begrijpt het niet. Wellicht zou jou wat onduidelijkheden kunnen geven (als je tijd hebt).

Klopt je berekening van de schwarzschildradius in je Python progje wel?

[OOOVincentOOO]

Goed gevonden!! Inderdaad een fout dat moet \(2\) zijn!!

Waarschijnlijk een typefout buiging hoek is \(\Delta \varphi=4MG/c^2R\) en \(r_s=2MG/c^2\). Stomme verwisseling \(4\) en \(2\)!

Ergens ben ik begonnen met \(R_s\) i.p.v. \(2MG/c^2\) te gebruiken. Jij bent het slachtoffer geworden!

Numeriek methode 1 formule (9):
$$\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} \,\,\,\,\,\,\, (9)$$
Correcte afbuig hoek en twee pieken:

Hier kan ik \(x^2/r^2\) niet elimineren zit afgeleide tussenin. Voor meer analyse zie openingspost: Lichtbuiging Analyse Twee Pieken (Mathpages). Vergelijkbare functie (mathpages) uitgesplits laat ook een en twee pieken zien.

Code: Selecteer alles

#Open pyplot in separate interactive window
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt5')

#https://www.mathpages.com/rr/s8-09/8-09.htm
#https://www.mathpages.com/rr/s6-03/6-03.htm

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams

rcParams['axes.titlepad'] = 20 

widths = [10,10]
heights = [10]

fig= plt.figure(figsize=(20,10))

gs=fig.add_gridspec(1,2,width_ratios=widths, height_ratios=heights)

ax1=fig.add_subplot(gs[0,0])
ax2=fig.add_subplot(gs[0,1])


#Physical constants
M=1.989e30
G=6.67408e-11
c=3e8
R=696340000

#schwarzschild Radius
Rs=2*M*G/c**2

#Function c(r)
def fdphidr1(r,x):
    
    #met x^/r^2:
    phi=(1/(1-Rs/r)-0.5*(-3*(x**2/r**2)+1)/((x**2/r**2) *Rs/r  +1 -Rs/r))*Rs*R/r**3
    
    #zonder x^2/r^2
    phi=(1/(1-Rs/r)-0.5*(-3+1)/(Rs/r  +1 -Rs/r))*Rs*R/r**3
  
    return phi

y=R
x=np.linspace(-10,10,10000)
x=x*R
r=np.sqrt(y**2+x**2)

#Angular Distribution
dphidr=fdphidr1(r,x)

ax1.plot(x/R,dphidr,color="red", linewidth=0.5, label=r"(9)")
ax1.set_xlabel(r'$x$',fontsize=15)
ax1.set_ylabel(r'$d \phi /dx$',fontsize=15)

#Integrated deflection agngle
deflection=np.cumsum(dphidr)*(x[10]-x[9])
deflection=np.degrees(deflection)*3600

ax2.plot(x/R,deflection,color="red", linewidth=0.5, label=r"(9)")
ax2.set_xlabel(r'$x$',fontsize=15)
ax2.set_ylabel(r'$\int \phi dx$',fontsize=15)


ax1.legend(loc="upper right")
ax2.legend(loc="upper right")

Numeriek Methode 2 op basis \(c(r)\):
Middels:
$$\frac{\partial c(r)}{\partial y}= \frac{d \varphi}{dx}$$
Formule Puntje:
$$\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} c(r)$$
$$c(r)=\frac{1}{2} \ln \left (\frac{(1 - \frac{r_s}{r})^2}{ \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} } \right )$$

Correcte afbuig hoek met \(x^/r^2\) (twee pieken):

Correcte afbuig hoek zonder \(x^/r^2\) (een piek):

Formule Mathpages:
$$c(r)=\sqrt{\frac{1-\frac{r_s}{r}}{1+\frac{r_s}{r}\frac{x^2}{r^2}\left(
\frac{1}{1-r_s/r} \right)}}$$
Correcte afbuig hoek met \(x^/r^2\) (twee pieken):

Correcte afbuig hoek zonder \(x^/r^2\) (een piek):

Voor meer analyse zie: Lichtbuiging Analyse Twee Pieken (Mathpages). Een en twee pieken bij restant: \(x^2/r^2\) door verwaarlozing \(dy\).

Code: Selecteer alles

#Open pyplot in separate interactive window
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt5')

#https://www.mathpages.com/rr/s8-09/8-09.htm
#https://www.mathpages.com/rr/s6-03/6-03.htm

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams

rcParams['axes.titlepad'] = 20 

widths = [5,5,5]
heights = [5,5]

fig= plt.figure(figsize=(30,10))

gs=fig.add_gridspec(2,3,width_ratios=widths, height_ratios=heights)

ax1=fig.add_subplot(gs[0,0])
ax2=fig.add_subplot(gs[0,1])
ax3=fig.add_subplot(gs[0,2])
ax4=fig.add_subplot(gs[1,0])
ax5=fig.add_subplot(gs[1,1])
ax6=fig.add_subplot(gs[1,2])

#Physical constants
M=1.989e30
G=6.67408e-11
c=3e8
R=696340000

#schwarzschild Radius
Rs=2*M*G/c**2

#Function c(r)
def cxy(x,y):
    
    #mathpages methode c(r)
    nomerator=1-2*M*G/((c**2)*np.sqrt(x**2+y**2))
    denomerator= 1+    (2*M*G/(c**2))*1/np.sqrt(x**2+y**2)   *  x**2/(x**2+y**2)  * 1/ ( 1-(2*M*G/(c**2))/np.sqrt(x**2+y**2))
    #cr=np.sqrt(nomerator/denomerator)
    
    #Approx Methode Puntje
    cr=0.5*np.log((1-Rs/np.sqrt(x**2+y**2))**2/( Rs/np.sqrt(x**2+y**2)+1 -Rs/np.sqrt(x**2+y**2)   ))
    
    return cr

size=5000

#Create 2D mesh and calculate c(r)
x=np.linspace(-30,30,size)
x=R*x

y=np.linspace(1,6,size)
y=R*y

X,Y =np.meshgrid(x,y)
z=cxy(X,Y)

#Function c(r) as function of x and y
cf1=ax1.contourf(x/R,y/R,z,levels=15, cmap='seismic',alpha=1)
fig.colorbar(cf1,ax=ax1, location="bottom")
ax1.set_title(r'$c(\sqrt{x^2+y^2})$',fontsize=20)
ax1.set_xlabel(r'$x$',fontsize=15)
ax1.set_ylabel(r'$y$',fontsize=15)

#Gradient for d/dy this is dTheta/dx
dcdy = np.gradient(z, axis=0)/(y[5]-y[4])
cf2=ax2.contourf(x/R,y/R,dcdy,levels=15, cmap='seismic',alpha=1)
fig.colorbar(cf2,ax=ax2, location="bottom")
ax2.set_title(r'$\partial c / \partial y = d \Theta /dx$',fontsize=20)
ax2.set_xlabel(r'$x$',fontsize=15)
ax2.set_ylabel(r'$y$',fontsize=15)

#Integrate Theta over z
intz=np.cumsum(dcdy,axis=1)*(x[5]-x[4])
cf3=ax3.contourf(x/R,y/R,np.degrees(intz)*3600,levels=15, cmap='seismic',alpha=1)
fig.colorbar(cf3,ax=ax3, location="bottom")
ax3.set_title(r'$\int \Theta \ dx$',fontsize=20)
ax3.set_xlabel(r'$x$',fontsize=15)
ax3.set_ylabel(r'$y$',fontsize=15)

#Function c(r) as function of x and y
for i in range(np.size(y)):
    
    if 5*(i/size)%1==0:
     
        cr=z[i]
        ax4.plot(x/R,cr,color="black", linewidth=0.5, label="x=" + str(x[i]))
        
        ax4.annotate("y=" + str(round(y[i]/R,2)), xy=( x[np.argmin(cr)]/R , np.min(cr)),ha='center', xycoords='data', xytext=(-100 , -5), textcoords='offset points',rotation = 0,arrowprops=dict(arrowstyle='-',lw=0.2))
        
        ax4.set_xlabel(r'$x$',fontsize=15)
        ax4.set_ylabel(r'$c(\sqrt{x^2+y^2})$',fontsize=15)
        
        ax4.ticklabel_format(axis='y', style='sci', scilimits=(0,0))

#Gradient for d/dy this is dTheta/dx        
p=0
for i in range(np.size(y)):

    if 5*(i/size)%1==0:
     
        
        dphidy=dcdy[i]
        ax5.plot(x/R,dphidy,color="black", linewidth=0.5, label="x=" + str(x[i]))

        ax5.annotate("y=" + str(round(y[i]/R,2)), xy=( x[np.argmax(dphidy)]/R , np.max(dphidy)),ha='center', xycoords='data', xytext=(((-1)**p)*100 ,-5), textcoords='offset points',rotation = 0,arrowprops=dict(arrowstyle='-',lw=0.2))
        
        ax5.set_xlabel(r'$x$',fontsize=15)
        ax5.set_ylabel(r'$d \Theta /dx$',fontsize=15)
        p+=1
        print(p)

#Integrate Theta over z       
for i in range(np.size(y)):
    
    if 5*(i/size)%1==0:
     
        phidx=np.degrees(intz[i])*3600
        ax6.plot(x/R,phidx,color="black", linewidth=0.5)
        
        ax6.annotate("y=" + str(round(y[i]/R,2)), xy=(  x[np.argmax(phidx)]/R , np.max(phidx)),ha='center', xycoords='data', xytext=(-50 ,-10), textcoords='offset points',rotation = 0,arrowprops=dict(arrowstyle='-',lw=0.2))
        
        ax6.set_xlabel(r'$x$',fontsize=15)
        ax6.set_ylabel(r'$\int \Theta \ dx$  [arcsec.]',fontsize=15)     

theta=4*G*M/(R*c**2)
theta=np.degrees(theta)*3600
ax6.plot([x[0]/R,x[-1]/R],[theta,theta],color="red", linewidth=0.5, label=r"$\dfrac{4GM}{Rc^2}$")
ax6.legend(loc="lower right")

Vraag:
Waarom geven beide \(c(r)\) dezelfde oplossingen voor \(\ln\) en \(\sqrt{}\)? Misschien probeer ik straks de afgeleide te nemen naar \(\partial / \partial y\) deze zouden dan gelijk zijn?

Puntje:
$$c(r)=\frac{1}{2} \ln \left (\frac{(1 - \frac{r_s}{r})^2}{ \frac{r_s}{r} \frac{x^2}{r^2} + 1 - \frac{r_s}{r} } \right )$$
Mathpages:
$$c(r)=\sqrt{\frac{1-\frac{r_s}{r}}{1+\frac{r_s}{r}\frac{x^2}{r^2}\left(
\frac{1}{1-r_s/r} \right)}}$$

Nu pauze deze post was HEEEEEEEL erg veel werk!!!!

[Professor Puntje]

Mooi!

Vermenigvuldig teller en noemer bij MathPages met 1 - r_s/r dan lijkt het er al meer op.

[Professor Puntje]

Mijn versie van Huygens' principe is preciezer dan die van MathPages:

Even herhalen wat MathPages met Huygens' principe doet. Dat komt op het volgende neer:

Uit bovenstaand plaatje van een infinitesimaal voortschrijdend golffront zien we dat:

\(\)

\( \mathrm{d} \varphi = \tan( \mathrm{d} \varphi ) \)

\(\)

\( \mathrm{d}\varphi = \frac{ (x_{R+dR}(t + \mathrm{d}t) - x_{R+dR}(t)) \, - \, (x_R(t + \mathrm{d}t) - x_R(t)) )}{\mathrm{d}R} \)

\(\)

\( \mathrm{d}\varphi = \frac{ \frac{\mathrm{d} x_{R+dR}}{\mathrm{d}t} \cdot \mathrm{d}t \, - \, \frac{\mathrm{d} x_R}{\mathrm{d}t} \cdot \mathrm{d}t }{\mathrm{d}R} \)

\(\)

\( \mathrm{d}\varphi = \frac{ \frac{\mathrm{d} x_{R+dR}}{\mathrm{d}t} \, - \, \frac{\mathrm{d} x_R}{\mathrm{d}t} }{\mathrm{d}R} \cdot \mathrm{d}t \)

\(\)

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

\(\)

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

\(\)

\( \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) \)

\(\)

Let wel dat hier sprake is van een benadering, en dat Huygens' principe zoals op MathPages toegepast nog in de race is als mogelijke (deel)oorzaak van de twee pieken.

Op MathPages mist de noemer van (7), dat is een verwaarlozing die ik niet toepas.

[Professor Puntje]

Verder heb ik voor het kwadraat van de horizontale lichtsnelheid:

\(\)

\( \frac{dx^2}{dt^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 \,\,\,\,\,\,\,\, (8) \)

[OOOVincentOOO]

Ik snap het nog niet allemaal:

Uit mathpages (twee pieken: Mathpages 8-09) en en uit jouw afleiding en deels Wiki (heb allen vergeleken, nagelopen en stemmen overeen) dan komt de volgende complete matrix uit Schwarzschild en vind dat eigenlijk de meest overzichtelijke methode. Deze staat ook in Mathpages 6-06 (waar 8-09 naar ref.):

$$g=\begin{bmatrix}\boxed{1} & 0 & 0 & 0 \\ 0 & \boxed{ -1} & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \end{bmatrix} - \frac{r_s}{r} \begin{bmatrix}\boxed{1} & 0 & 0 & 0 \\ 0 & \boxed{\kappa x^2} & \kappa xy & \kappa xz \\ 0 & \kappa yx & \kappa y^2 & \kappa yz \\ 0 & \kappa zx & \kappa zy & \kappa z^2 \end{bmatrix} \\ \kappa=\frac{1}{r^2 \left(1- \frac{r_s}{r} \right)} $$
Hoe ik begrijp, men maakt alleen gebruik van de termen met de box. Dus zonder \(dy\) met als resultaat een niet symmetrisch model door: \(x^2/r^2\) omdat de rest allemaal als functie van \(r\) is. Door: \(x^2/r^2\) ontstaan de pieken (zie mijn eerdere reply) hoe ik begrijp en kan zelfs \(x^2/r^2=1\) gesteld worden en krijg je maar een piek (empirisch gevonden). Met en zonder \(x^2/r^2\) dezelfde afbuighoek.
$$g_{tt}=\left( \frac{d \tau}{dt} \right)^{2}=1-\frac{r_s}{r}
\\
g_{xx}=\left( \frac{d \tau}{dx} \right)^{2}=-1-\frac{r_s}{r}\frac{x^2}{r^2}\left(
\frac{1}{1-r_s/r} \right) $$
$$c(r)=\frac{dx}{dt}=\sqrt{-\frac{g_{tt}}{g_{xx}}}
\\
c(r)=\sqrt{\frac{1-\frac{r_s}{r}}{1+\frac{r_s}{r}\frac{x^2}{r^2}\left(
\frac{1}{1-r_s/r} \right)}}$$
En tenslotte de hoekverdeling \(d \varphi /dx\) en \(r^2=x^2+y^2\) waar het om gaat (welke ik numeriek oplos):
$$\frac{\partial c(r)}{\partial y}= \frac{d \varphi}{dx}$$

Bovenstaand is hoe ik het begrijp met \(g_{tt}\) en \(g_{xx}\) de afgeleiden.

In welke van bovenstaande stappen word dan een benadering gemaakt volgens jouw?
Zou je dat bovenstaand kunnen aangeven?
Dan kan ik wat verder studeren!

[Professor Puntje]

Dit maakt het onnodig ingewikkeld. Ik ga uit van de schwarzschildmetriek zoals die normalerwijze wordt uitgeschreven, en transformeer dan naar het xy-frame. Aangezien het licht in één vlak beweegt waarin ook het middelpunt van de zon ligt kun je dan de z-as daar haaks op kiezen en daarmee verdwijnen dan alle termen met z en dz. Dat is nog geen benadering maar enkel een handige coördinatenkeuze. Vervolgens wordt dy = 0 gesteld wat wel een benadering is. En dat geeft dan mijn formule:

\(\)

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

[Professor Puntje]

Onderstaande is een benadering:

\(\)

\( \frac{\partial c(r)}{\partial y}= \frac{d \varphi}{dx} \)

\(\)

Daar heb ik een preciezere formule voor gebruikt, zie mijn eerdere antwoord.