Python en Jacobi

[CoenCo]

@Vincent:
Dat is een setting in spyder:

Preferences -> IPython console -> Graphics -> Graphics Backend -> zet op "Automatic" (staat nu waarschijnlijk op "Inline")

[OOOVincentOOO]

Hallo,

Momenteel maintenance op werk. Ben bezig met Python training .

Hallo Puntje,
Met scipi: .ellipk kan je de startconditie berekenen zonder gebruik van wolfram [Scipy]. Hier komt: 1.5707996366183763. Ik weet niet precies wat ik allemaal doe maar het schijnt te werken.

Code: Selecteer alles

import scipy.special as sps

#Solve angle start conditions
h=8.428356322064124e-06
angle=sps.ellipk(h)
print('angle: ' + str(angle))

Onderstaand heb ik nog wat gespeeld en een slider toegevoegd voor hoek \(\alpha\) de zoom factor. Ik weet niet of de code bij jouw gaat werken. Als het goed is opent een apart window met interactieve grafiek.

Ik begrijp iets niet. Indien de hoek range zeer klein is krijg ik een halve circel. Misschien begrijp jij wat er gebeurt? Zijn de getallen dan te klein? [EDIT: Zie zojuist dat benaming Phi verkeerd was lijkt nu te werken, code gupdate]

Betreffende de code. Ik heb het zo simpel mogelijk proberen te houden. Her en der iets van commentaar. Hopelijk krijg je zin om the pimpen (steal the best and make the rest).

Code: Selecteer alles

#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 matplotlib.widgets import Slider
fig, ax1= plt.subplots(figsize=(9, 5))

# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.2, -0.01,0.5, 0.05], facecolor=axcolor)
s_alpha = Slider(axalpha, 'alpha: ', 0, 1, valinit=1)


#Set constants note G is divided by c^2    
M=1.989e30
G=6.67408e-11/((3e8)**2)
Ro=7e8   

#Calculate: e1, e2, e3
esqrt=np.sqrt((Ro-2*M*G)*(Ro+6*M*G))
e1=(Ro-2*M*G+esqrt)/(4*M*G*Ro)
print('e1: ' + str(e1))

e2=1/Ro
print('e2: ' + str(e2))

e3=(Ro-2*M*G-esqrt)/(4*M*G*Ro)
print('e3: ' + str(e3))

Tau=np.sqrt(M*G*(e1-e3)/2)
print('Tau: ' + str(Tau))

h=(e2-e3)/(e1-e3)
print('h: ' + str(h))

#Solve angle start conditions
angle=sps.ellipk(h)
print('angle: ' + str(angle))

#Calculate Sigma Start Condition
Sigma=-Tau*np.pi/2-angle
print('sigma: ' + str(Sigma))

#Set phi range default plot
phi=np.linspace(np.pi/2-1*np.pi/2,np.pi/2+1*np.pi/2,1000)

#Calculate Radius
v=Tau*phi+Sigma
sn,_,_,_=sps.ellipj(v,h)
r=1/(e3+(e2-e3)*(sn**2))

x=r*np.cos(phi)
y=r*np.sin(phi)

deflection, = ax1.plot(x,y)
ax1.grid()
  
#The following function updateplot runs whenever slider changes value
def updateplot(val):
    alpha = s_alpha.val
    phi=np.linspace(np.pi/2-alpha*np.pi/2,np.pi/2+alpha*np.pi/2,1000)

    v=Tau*phi+Sigma
    sn,_,_,_=sps.ellipj(v,h)
    r=1/(e3+(e2-e3)*(sn**2))
    
    x=r*np.cos(phi)
    y=r*np.sin(phi)
    ax1.clear()
    deflection, = ax1.plot(x,y)
    
    ax1.grid()   
    fig.canvas.draw_idle()

#The following code checkes if slider changes. This line is looped automatic by pyplot
s_alpha.on_changed(updateplot);

[OOOVincentOOO]

Hallo,

Nog een grapje toegevoegd. Ik heb een circel met radius zon toegevoegd. Voor normale setting is de deflektie niet zichtbaar. Dus x en y dezelfde schalering.

Ik heb de zon massa 25000 maal zo zwaar gemaakt (even vergeten of het physisch mogelijk is). Dan krijgt men een indruk van de stralengang tov. radius.

Iets minder programmeren maar meer analyseren.

Code: Selecteer alles

#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 matplotlib.widgets import Slider
fig, ax1= plt.subplots(figsize=(9, 5))

# Draw slider alpha
axcolor = 'lightgoldenrodyellow'
axalpha = plt.axes([0.2, -0.01,0.5, 0.05], facecolor=axcolor)
s_alpha = Slider(axalpha, 'alpha: ', 0, 1, valinit=1)


#Set constants note G is divided by c^2    
M=25000*1.989e30
G=6.67408e-11/((3e8)**2)
Ro=7e8   

#Calculate: e1, e2, e3
esqrt=np.sqrt((Ro-2*M*G)*(Ro+6*M*G))
e1=(Ro-2*M*G+esqrt)/(4*M*G*Ro)
print('e1: ' + str(e1))

e2=1/Ro
print('e2: ' + str(e2))

e3=(Ro-2*M*G-esqrt)/(4*M*G*Ro)
print('e3: ' + str(e3))

Tau=np.sqrt(M*G*(e1-e3)/2)
print('Tau: ' + str(Tau))

h=(e2-e3)/(e1-e3)
print('h: ' + str(h))

#Solve angle start conditions
angle=sps.ellipk(h)
print('angle: ' + str(angle))

#Calculate Sigma Start Condition
Sigma=-Tau*np.pi/2-angle
print('sigma: ' + str(Sigma))

#Set phi range default plot
phi=np.linspace(np.pi/2-1*np.pi/2,np.pi/2+1*np.pi/2,1000)

#Calculate Radius
v=Tau*phi+Sigma
sn,_,_,_=sps.ellipj(v,h)
r=1/(e3+(e2-e3)*(sn**2))

x=r*np.cos(phi)
y=r*np.sin(phi)

#Plot circle radius for sun
circle=plt.Circle((0,np.max(y)-Ro),Ro, color='orange')
ax1.add_patch(circle)

#Set axis equal so sun is round
ax1.axis('equal')
    
deflection, = ax1.plot(x,y)
ax1.grid()
  
#The following function updateplot runs whenever slider changes value
def updateplot(val):
    alpha = s_alpha.val
    phin=np.linspace(np.pi/2-alpha*np.pi/2,np.pi/2+alpha*np.pi/2,1000)

    v=Tau*phin+Sigma
    sn,_,_,_=sps.ellipj(v,h)
    r=1/(e3+(e2-e3)*(sn**2))
    
    x=r*np.cos(phin)
    y=r*np.sin(phin)
    ax1.clear()
    deflection, = ax1.plot(x,y)
    
    #Plot circle radius for sun
    circle=plt.Circle((0,np.max(y)-Ro),Ro, color='orange')
    ax1.add_patch(circle)
    
    #Set axis equal so Sun is round
    ax1.axis('equal')
    
    ax1.grid()   
    fig.canvas.draw_idle()

#The following code checkes if slider changes. This line is looped automatic by pyplot
s_alpha.on_changed(updateplot);

[Professor Puntje]

De code werkt bij mij ook.

Eerst eens even kijken wat je met sps.ellipk doet, en of ik dat in mijn eigen progje kan overnemen. Dat zou het weer wat eleganter maken als ik sigma niet zelf handmatig hoef te benaderen.

[Professor Puntje]

Hallo Puntje,
Met scipi: .ellipk kan je de startconditie berekenen zonder gebruik van wolfram [Scipy]. Hier komt: 1.5707996366183763. Ik weet niet precies wat ik allemaal doe maar het schijnt te werken.

Code: Selecteer alles

import scipy.special as sps

#Solve angle start conditions
h=8.428356322064124e-06
angle=sps.ellipk(h)
print('angle: ' + str(angle))

Dat klopt. Zie:

En:

Weer even terug naar:

\(\)

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

\(\)

Het is goed mogelijk dat we h² in plaats van h moeten nemen. Immers heeft men het bij ellipj over de parameter m in plaats van k. Daarom bepaal ik nu σ uitgaande van de aangepaste formule (6*):

\(\)

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

\(\)

Nu bekijken we de lichtbaan waarvoor r minimaal (d.w.z. r₀) is voor φ = π/2. Dat geeft:

\(\)

\( r_0 = \frac{1}{ e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h^2)} \)

\(\)

En wegens (2) dat:

\(\)

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

\(\)

\( e_2 = e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h^2) \)

\(\)

\( e_2 - e_3 = (e_2 - e_3) \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h^2) \)

\(\)

\( \mathrm{sn}^2( \tau \frac{\pi}{2} + \sigma ; h^2) = 1 \)

\(\)

Aangezien sn ook voor reële argumenten periodiek is zijn daar oneindig veel oplossingen voor. We nemen een simpel geval en bekijken:

\(\)

\( \mathrm{sn}( \tau \frac{\pi}{2} + \sigma ; h^2) = 1 \)

\(\)

Proberenderwijs met Spyder en ellipj kom ik op:

\(\)

\( \sigma = -\tau \frac{\pi}{2} + 1,57079965 \)

\(\)

Die laatste stap moet dus ook via sps.ellipk en m=k² kunnen. Ga even mijn progje aanpassen.

[Professor Puntje]

Hier het aangepaste progje:

Code: Selecteer alles

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

alpha = 8.42e-6

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(-alpha/2,np.pi+alpha/2,1000)

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

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

plt.plot(x,y,c=(1,0.2,0.5),lw=1)
plt.title('lichtbaan')
plt.axis('on')
plt.show()

[Professor Puntje]

Ik ga nu nog even bekijken of alpha ook automatisch berekend kan worden. Dat is wel zo fraai.

[Professor Puntje]

Uit (6*) zien we dat de asymptoten optreden voor:

\(\)

\( e_3 + (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h^2) = 0 \)

\(\)

\( (e_2 - e_3) \mathrm{sn}^2( \tau \varphi + \sigma ; h^2) = -e_3 \)

\(\)

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

\(\)

Zoiets is volgens WolframAlpha oplosbaar met de inverse van sn. Zie:

Maar of die inverse functie ook voor Python voorradig is?

[OOOVincentOOO]

Nice, even voor de duidelijkheid. Zonder jouw inzicht over plot range: 0 tot pi en de methode om sigma te bepalen was ik niet uitgekomen.

Ik ben een beetje aan het stoeien met numerieke integratie. Niet het oppervlakte maar de integraal functie

Methode 1:
Via numpy cumulatieve sommatie (deze heb ik eerder wel gebruikt):

Code: Selecteer alles

import numpy as np
a=[0,1,2,3,4,5]
ai=np.cumsum(a)
print(ai)

Methode 2:
Via Scipy

Code: Selecteer alles

import scipy.integrate as integrate
a=[0,1,2,3,4,5]
ai=integrate.cumtrapz(a)
print(ai)

x=[2,4,8,16,32]
y=[1,1,1,1,1]
ai=integrate.cumtrapz(y,x)
print(ai)

Ik heb moeilijkheden te begrijpen, de x-as is niet linear verdeelt. De scipy geeft de mogelijkheid ook een x waarde in te voeren. Bij numpy moet je dat "handmatig" doen. Volgens mij werkt scipy prima indien de x-as ingevoerd

Ik stoei wat verder met de integratie. Indien je via scipy integreert verlies je 1 element in het array. Ik kan je wat trukjes laten zien hoe hiermee om te. Mits je het wil weten.

[Professor Puntje]

Heb nu ook de "separate interactive window" toegevoegd.

[OOOVincentOOO]

Uit (6*) zien we dat de asymptoten optreden voor:

\(\)

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

\(\)

Zoiets is volgens WolframAlpha oplosbaar met de inverse van sn. Zie:

Maar of die inverse functie ook voor Python voorradig is?

Weet niet precies wat je aan het koken bent. Maar is dit niet hetzelfde als:

Code: Selecteer alles

import scipy.special as sps

#Solve angle start conditions
h=8.428356322064124e-06
angle=sps.ellipk(h)
print('angle: ' + str(angle))

Of zoek je de oplossing met complexe getallen?

[Professor Puntje]

Dat codefragment heb ik net al in mijn aangepaste progje toegepast. Maar dat geeft nog geen antwoord op de vraag wat de totale lichtafbuiging is. Hoe ik het begrijp is het effect van sigma dat door veranderingen in de waarde van sigma de hele grafiek rond de oorsprong roteert. Dat willen we niet, want we hebben graag dat de lichtbaan de gebruikelijk stand heeft waarbij het licht voor φ=π/2 de zon bovenlangs (met r = r₀) het dichtst benadert. Die eis levert de waarde van sigma op. Maar daarmee weten we de totale lichtafbuiging nog niet. Die totale lichtafbuiging heb ik al wel gevonden door de geplotte range van φ stapsgewijze te vergroten tot op het moment dat je onzinnige plots krijgt. Maar dat is een kunstgreep die beter nog door een nette berekening kan worden vervangen.

[CoenCo]

Nice, even voor de duidelijkheid. Zonder jouw inzicht over plot range: 0 tot pi en de methode om sigma te bepalen was ik niet uitgekomen.

Ik ben een beetje aan het stoeien met numerieke integratie. Niet het oppervlakte maar de integraal functie

Methode 1:
Via numpy cumulatieve sommatie (deze heb ik eerder wel gebruikt):

Code: Selecteer alles

import numpy as np
a=[0,1,2,3,4,5]
ai=np.cumsum(a)
print(ai)

Methode 2:
Via Scipy

Code: Selecteer alles

import scipy.integrate as integrate
a=[0,1,2,3,4,5]
ai=integrate.cumtrapz(a)
print(ai)

x=[2,4,8,16,32]
y=[1,1,1,1,1]
ai=integrate.cumtrapz(y,x)
print(ai)

Ik heb moeilijkheden te begrijpen, de x-as is niet linear verdeelt. De scipy geeft de mogelijkheid ook een x waarde in te voeren. Bij numpy moet je dat "handmatig" doen. Volgens mij werkt scipy prima indien de x-as ingevoerd

Ik stoei wat verder met de integratie. Indien je via scipy integreert verlies je 1 element in het array. Ik kan je wat trukjes laten zien hoe hiermee om te. Mits je het wil weten.

Code: Selecteer alles

scipy.integrate.cumtrapz(y,x, initial=0)

Levert hetzelfde aantal elementen op als je er in stopt. De initial is dan de startwaarde.

Waar zit je precies mee verder?
Cumsum is precies wat het zegt:de cumulatieve som van alle elementen, waarbij ook alle tussenresultaten gegeven worden.

Cumtrapz is de cumulutatieve som van alle trapeziums tussen de gegeven punten. Voor het oppervlak heb je natuurlijk de breedte (x) en hoogte (y1, y2) nodig. Daarnaast passen er tussen 10 elementen, maar 9 trapeziums, daarom krijg je er standaard eentje minder terug.
Zie ook trapezoidal integration: https://en.m.wikipedia.org/wiki/Trapezoidal_rule

[OOOVincentOOO]

@Coenco,
Dat was voor puntje dat hij niet hoeft ze zoeken naar mogelijke methoden. De methoden zijn mij al veel eerder bekend. Betreffende topic: de x-as is niet verdeeld in gelijke bins want de hoek wordt omgerekend naar de x-waarde. Daar ben ik wat mee aan het puzzelen.

@Puntje,
Ik begrijp het nog niet helemaal. Mijn versie van vanmiddag had een bereik voor alpha van: 0 tot pi (en met de slidebar kan je inzoomen tot bijvoorbeeld jouw eerste variant plot). Wil jij een model wat ook hoeken beschrijft groter dan pi?

[Professor Puntje]

Heb je interactieve plot nog weer even bekeken, en volgens mij bekijk je φ van 0 tot π (in het geval 100%) en voor kleinere percentages een kleiner interval maar wel steeds zo dat φ = π/2 in het midden zit. Dat zit dan allemaal boven of op de x-as. Maar in je plots zie ik ook negatieve y-waarden. Dus ik begrijp niet wat daar gebeurt.