Python en Jacobi

[Professor Puntje]

Ja - werkt bij mij en geeft mooie plaatjes. Bij sommige instellingen krijg je wel vreemde curves. Maar als de straal ten opzichte van de massa te klein wordt krijg je een zwart gat. Dat kun je nagaan door r₀ met r_s met te vergelijken.

[Marko]

Graag wat specifieker. Aan vage suggesties hebben we niets.

Concreet: In veel gevallen kan een eerste orde benadering die tot een simpele analytische oplossing tot een nauwkeurigere benadering van de werkelijkheid komen dan een exacte analytische oplossing die voor het daadwerkelijke uitrekenen tig numerieke rekenstappen vereist.

Strikt genomen moet je het natuurlijk sowieso omdraaien en aantonen dat de benaderingen in de gebruikte methode acceptabel zijn en niet tot artefacten kunnen leiden, voor je een conclusie trekt zoals je die trekt. Zaken waar je op moet letten:

1. Die elliptische functies zijn kennelijk een analytische oplossing van een DV. Maar de waarde van die oplossing wordt binnen die Python functie nog steeds numeriek berekend en je zou uit de documentatie moeten halen met welke nauwkeurigheid dat is en onder welke voorwaarden

2. Die vraag over nauwkeurigheid en voorwaarden geldt in feite voor allerlei numerieke bewerkingen die in deze aanpak worden gebruikt. Dus bij het bepalen van de wortel, bij differentiëren, enzovoort.

[OOOVincentOOO]

Ja - werkt bij mij en geeft mooie plaatjes. Bij sommige instellingen krijg je wel vreemde curves. Maar als de straal ten opzichte van de massa te klein wordt krijg je een zwart gat. Dat kun je nagaan door r₀ met r_s met te vergelijken.

Inderdaad er is nog een gedeeld door nul situatie bij de radius. Ik heb die criteria eerst in code gehad maar eruit gehaald voor overzicht.

Datzelfde radius issue had ik ook met genormaliseerde grootheden met R=1, G=1 enzo. Dan kwam ik niet op de factor 4 uit welke in de deflectie lichtbuiging formule staat.

[Professor Puntje]

@ Marko

Ja - dat klopt. De benaderingen die hier nog inzitten komen uit de toepassing van numerieke methoden (oftewel het gebruik van Python). Ik had de zaak ook liever helemaal analytisch aangepakt, maar dat is me niet gelukt. De formules werden daar te ingewikkeld voor. Hoe nauwkeurig is Python? Wellicht kunnen anderen met meer ervaring met dat programma daar beter op antwoorden. Uit de berichten hier op het Wetenschapsforum begreep ik dat het nauwkeurig genoeg is voor een simulatie als in dit topic.

[OOOVincentOOO]

Volgens Theorie [Wiki]:
$$\theta ={\frac {4GM}{rc^{2}}}$$
$$\theta ={\frac {(4)(6.67408e-11)(1.989e30)}{(69634e4)(3e8)^{2}}}=1.747621839''$$
$$\frac{\theta}{2} = 0,873810919''$$

Volgens Python met Jacobi Elliptic (maximale precisie 100% theoretische bereik, kan geoptimaliseerd worden met meer datapunten):
$$\frac{\theta}{2} = 0,873815''$$

100% Bereik:

Voor aardse begrippen hetzelfde zou ik zeggen velen malen kleiner dan de meetfout.

Indien men inzoomt dus data weggooit bij Jacobi wordt de fout groter. Men neemt immers niet de limiet tot oneindig zoals: mathpages en Eddington.

50% Bereik:

[Professor Puntje]

In principe zouden we hier zelfs een exacter antwoord dan met de bekende formule \( \theta ={\frac {4GM}{rc^{2}}} \) moeten kunnen krijgen door de nauwkeurigheid van de numerieke methoden voldoende op te voeren. Die genoemde formule is immers maar een benadering. Maar hoe exact ons resultaat feitelijk is hangt ervan af hoe precies Python (en onze computers) rekenen. En dat laatste weet ik niet.

[OOOVincentOOO]

Edit:
Voor radius zon heb ik zowel in Python als theoretische formule: 69634e4 [m] gebruikt. In de code welke ik geplaats heb staat een afgeronde waarde. Dus opletten welke constanten men invult.

Reactie:
Op dit moment is onze Phi schaal linear verdeelt maar de omgerekende x schaal niet. Voor de uiterste punten waar de deflectie berekend wordt is de afwijking het grootste daar zijn de minste datapunten. Men zou een Phi (volgens mij arctan) schaal verdeling kunnen nemen waarbij er meer datapunten zijn in de uiterste punten. En op de tweede plaats inderdaad de precisie van aantallen decimalen. Echter allemaal verwaarloosbaar met de meetnauwkeurigheid in de menseljke wereld.

[Professor Puntje]

Maar die uiterste punten zijn voor de pieken-kwestie niet zo interessant. Daarvoor zal het - vermoed ik - eerder gaan om het aantal decimalen en de grootte van de stapjes bij het differentiëren. En om de nauwkeurigheid van de gebruikte functies.

[OOOVincentOOO]

Ohh, ik dacht dat het ging op de nauwkeurigheid Jacobi aanpak van de deflectie hoek t.o.v. "theorie". De stap grote differentiëren is precies waar ik het over heb. In de uiterste waarden van x zijn momenteel de stapjes dx het grootste door niet lineare verdeling (veroorzaakt door lineare schaalverdeling phi omrekend vervormen de x schaalafstanden x=r.cos(phi) dx).

[OOOVincentOOO]

Maar die uiterste punten zijn voor de pieken-kwestie niet zo interessant. Daarvoor zal het - vermoed ik - eerder gaan om het aantal decimalen en de grootte van de stapjes bij het differentiëren. En om de nauwkeurigheid van de gebruikte functies.

Ik heb net een simulatie gedraaid met slecht 5 punten voor phi. Ook hier: \(\theta/2=0.873815''\). Dus mijn eerste intuitie betreffende dx op uitersten (niet lineariteit) lijkt verwaarloosbaar.

Edit:
Aanpak Jacobi Elliptic is natuurlijk geen bewijs dat de deflektie een limiet is (convergeerd)!

[OOOVincentOOO]

Nu moet ik stoppen. Ik zie net dat er een spelfout in staat: Jabobi ipv Jacobi . Daarom wil ik even eerherstel doen:

Versie 11.0: Spelling Jacobi/Distribution en massa zon preciezer.

Code: Selecteer alles

#Version 11.0 Schwarzschild, Jacobi Elliptic
#https://www.wetenschapsforum.nl/viewtopic.php?f=85&t=212427

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

import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sps
from scipy import optimize
from matplotlib.widgets import Slider, CheckButtons                    

fig= plt.figure(figsize=(15, 6))  

ax1 = plt.subplot2grid((3, 3), (0, 0),rowspan=2)
ax2 = plt.subplot2grid((3, 3), (0, 1),rowspan=2)
ax3 = plt.subplot2grid((3, 3), (0, 2),rowspan=2)
ax4 = plt.subplot2grid((3, 3), (2, 0),colspan=3)

ax4.axis('off') 

plt.suptitle('Schwarzschild, Jacobi Elliptic', fontsize=22) 
                 
# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.25, 0.05,0.4, 0.04], facecolor=axcolor)
s_zoom = Slider(axalpha, 'zoom: ', 0, 1, valinit=1)

axsm = plt.axes([0.25, 0.1,0.4, 0.04], facecolor=axcolor)
s_sm = Slider(axsm, 'solar masses: ', 1,100000, valinit=1)

radius = plt.axes([0.25, 0.15,0.4, 0.04], facecolor=axcolor)
s_sr = Slider(radius, 'radius: ', 0.1, 5, valinit=1)

rax = plt.axes([0.75, 0.05, 0.15, 0.15])
check = CheckButtons(rax, ('xy equal', 'sun'), (False, False))

#Set constants note G is divided by c^2    
global M
global G
global Ro

M=1.989e30
G=6.67408e-11/((3e8)**2)
Ro=696340000

def e1(sm,sr):
    esqrt=np.sqrt((sr*Ro-sm*2*M*G)*(sr*Ro+sm*6*M*G))
    e1=(sr*Ro-sm*2*M*G+esqrt)/(sm*4*M*G*sr*Ro)
    return e1

def e2(sr):
    e2=1/(sr*Ro)
    return e2

def e3(sm,sr):
    esqrt=np.sqrt((sr*Ro-sm*2*M*G)*(sr*Ro+sm*6*M*G))
    e3=(sr*Ro-sm*2*M*G-esqrt)/(sm*4*M*G*sr*Ro)
    return e3   
    
def r(phi):
    sm = s_sm.val
    sr = s_sr.val
     
    e1v=e1(sm,sr)
    e2v=e2(sr)
    e3v=e3(sm,sr)
        
    tau=np.sqrt(sm*M*G*(e1v-e3v)/2)
    h = (e2v-e3v)/(e1v-e3v)
    sigma = -tau*np.pi/2 - sps.ellipk(h)   
    sn,_,_,_=sps.ellipj(tau*phi+sigma, h)
    r=e3v+(e2v-e3v)*sn**2
    return r

#The following function updateplot runs whenever slider changes value
def updateplot(val):
    

    clearsetplot()

    #Solar masses/radius and magnification alpha
    zoom = s_zoom.val
    sm = s_sm.val
    sr = s_sr.val

    #Find root angle for non-solutions
    root = optimize.newton(r,0)
    alphamax =np.abs(root) # beveiliging tegen singulariteiten

    #Set angle range and calculate radius
    phi= np.linspace(alphamax/2,np.pi-alphamax/2,500)
    phi=phi[int(0.5*(1-zoom)*500):int(0.5*(1+zoom)*500)]
   
    rv=1/r(phi)
    
    x=rv*np.cos(phi)
    y=rv*np.sin(phi)
    
    #plot lightpath curve
    ax1.plot(x,y)
    ax1.grid()
            
    #Angular Deflection
    dx=np.diff(x)
    dy=np.diff(y)
    dydx=dy/dx
    da=np.arctan(dydx)
    dam=np.max(da)
    ax2.plot(x[:-1],-da)
    ax2.grid()
    
    #Angular Distribution
    ddy=np.diff(da)
    ddydx=ddy/dx[:-1]
    ax3.plot(x[:-2],-ddydx)        
    ax3.grid()       
     
    text=(  str('{:10.6f}'.format(dam))      +   ' $[rad]$\n' +
            str('{:10.6f}'.format(np.degrees(dam)))      +   '$^{\circ}$\n' +
            str('{:10.6f}'.format(np.degrees(dam)*3600)) +  "$''$")
    ax1.text(0.01, 0.75,text, transform=ax1.transAxes, fontsize=14)
    
    #scaling method and draw sun
    plotoptions(x,y,sr)
    fig.canvas.draw_idle()
    
#This function set scaling axis and drawes sun
def plotoptions(x,y,sr):
    i=0
    for r in check.get_status():
       
        if i==0:
            #Set scaling auto or equal
            if r==False:
                ax1.axis('auto')
                range=np.max(y)-np.min(y)
                ax1.set_ylim([np.min(y)-range*2/100,np.max(y)+range*2/100])
                range=np.max(x)-np.min(x)
                ax1.set_xlim([-np.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax2.set_xlim([-np.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax3.set_xlim([-np.max(x)-range*2/100,np.max(x)]+range*2/100)
                
                
            else:
                ax1.axis('equal')
                range=np.max(x)-np.min(x)
                ax1.set_xlim([-np.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax2.set_xlim([-np.max(x)-range*2/100,np.max(x)]+range*2/100)
                ax3.set_xlim([-np.max(x)-range*2/100,np.max(x)]+range*2/100)
                                
        if i==1:
            #Draw Sun yes/no
              if r==True:
                circle=plt.Circle((0,np.max(y)-sr*Ro),Ro, color='orange',label='Label')
                ax1.add_patch(circle)
                
                circle=plt.Circle((0,np.max(y)-sr*Ro),sr*Ro,fill=False, ls='--',ec='black')
                ax1.add_patch(circle)
                
        i=i+1

def clearsetplot():
    ax1.clear()
    ax1.set_title('Light Path', fontsize=14)
    ax1.set_xlabel('x [m]', fontsize=14)
    ax1.set_ylabel('y [m]', fontsize=14)
    ax1.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
    
    ax2.clear()
    ax2.set_title('Angular Deflection', fontsize=14)
    ax2.set_xlabel('x [m]', fontsize=14)
    ax2.set_ylabel('deflection angle  [rad]', fontsize=14)
    ax2.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
    
    ax3.clear()
    ax3.set_title('Angular Distribution', fontsize=14)
    ax3.set_xlabel('x [m]', fontsize=14)
    ax3.set_ylabel('angular change [rad/m]', fontsize=14)
    ax3.ticklabel_format(axis='both', style='sci', scilimits=(0,0))
       
#Display plot n startup.
updateplot(1)

#The following code checkes if slider changes. This line is looped automatic by pyplot
s_sm.on_changed(updateplot);
s_sr.on_changed(updateplot);
s_zoom.on_changed(updateplot);
check.on_clicked(updateplot);

[Professor Puntje]

De totale afbuiging is niet controversieel. Zowel MathPages, als Eddington, als wijzelf komen op (vrijwel) dezelfde waarde uit. De discussie was over de vraag hoeveel pieken er voorkomen in de momentane afbuigingstoename dφ/dx als functie van x. MathPages vond er 2, bij Eddington zijn er 3 (of meer), en bij ons nu 1. Je kunt met een verschillend aantal pieken dus toch de juiste totale afbuiging vinden. Zowel MathPages als Eddington gebruiken benaderingen. Bij ons zitten de enige benaderingen in het rekenwerk van Python. Het lijkt mij heel onaannemelijk dat Python ons door onnauwkeurigheden van het rekenwerk een onjuist aantal pieken voorschotelt. Maar een hard bewijs voor de nauwkeurigheid van Python in ons specifieke geval heb ik niet. Daarvoor weet ik er als beginner met Python te weinig vanaf.

[OOOVincentOOO]

Ook hebben we geen bewijs dat de deflectie convergeerd met de Jacobi methode. Of weet jij een handigheid?

[HansH]

Ik vind de code lastig te lezen omdat het geen nette formules zijn zoals in mathcad.
Als ik begrijp welke formule ik moet gebruiken en wat je numeriek moet oplossen kan ik het wel eens in mathcad stoppen en dan kunnen we het vergelijken met python

[Professor Puntje]

Dit kan met de schwarzschildstraal \( r_s \) van de zon vereenvoudigd worden tot:

\( e_1 = \frac{r_0 - r_s + \sqrt{(r_0 - r_s)(r_0 + 3r_s)}}{2 r_s r_0} \,\,\,\,\,\,\, (1) \)

\(\)

\( e_2 = \frac{1}{r_0} \,\,\,\,\,\,\, (2) \)

\(\)

\( e_3 = \frac{r_0 - r_s - \sqrt{(r_0 - r_s)(r_0 + 3r_s)}}{2 r_s r_0} \,\,\,\,\,\,\, (3) \)

\(\)

\( \tau = \sqrt{ \frac{ r_s (e_1 - e_3)}{ 4 }} \,\,\,\,\,\,\, (4) \)

\(\)

\( h = \sqrt{ \frac{ e_2 - e_3}{ e_1 - e_3 }} \,\,\,\,\,\,\, (5) \)

\(\)

\( r = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h)} \,\,\,\,\,\,\, (6) \)

\(\)

Verder is \( r_0 \) de kleinste afstand r van de lichtbaan tot het centrum van de zon in de gebruikelijke schwarzschildcoördinaten.

Dit zijn de basisformules, alleen moest (6) later aangepast worden vanwege notatie-kwesties voor de sn-functie. Dat zul je voor jouw programma moeten uitzoeken welke notatie daarin voor de sn-functie gebruikt wordt.