Valversnelling in Chimay

[HansH]

als ik dat een maal heb dan kan mathcad dat heel snel uitrekenen en vervolgens kun je dan de volgende stappen toevoegen:
-dichtheidsfunctie toevoegen
-zwaartekrachtsvector bepalen in punt p(x,y,z) voor elk massa elementje in de integraal
-alle zwaartekrachts vectoren samen integreren tot 1 complete vector

[HansH]

ik ben er al uit:

3dintegraal1 460 keer bekeken

dus kan nu de volgende stappen doen. volume berekenen duurt maar ca 0.5 sec met convergentietolerantie op 10^-5, dus kan heel nauwkeurig rekenen.

[Xilvo]

als ik dat een maal heb dan kan mathcad dat heel snel uitrekenen

Of het sneller is dan mijn methode betwijfel ik. Maar het hangt van veel factoren af, zoals het aantal massapunten dat je gebruikt. Onder de kap moeten toch dezelfde berekeningen worden uitgevoerd, als je het numeriek doet.

-dichtheidsfunctie toevoegen

Als je weet hoe dat moet hoef je je niet over de integratiegrenzen te bekommeren, dan zet je de dichtheid op nul voor punten buiten de ellipsoïde.

[wnvl1]

Volume berekenen gaat in scipy snel. Nu probeer ik de gravitatie potentiaal te berekenen.

Code: Selecteer alles

import numpy as np
from scipy.integrate import tplquad

a=6356000 #rstraal bij pool
b=6378000 #rstraal bij evenaar
c=6378000 #rstraal bij evenaar

G=6.6743e-11 

r_loc=1 #positie Chimay
phi_loc=0 #positie Chimay
theta_loc=np.radians(40) #positie Chimay

rho = 5515 #massa dichtheid aarde

tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/np.absolute((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc) )**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc)  )**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc)  ) ),
                0, 1,
                lambda theta: 0, np.pi,
                lambda theta, phi: 0, 2*np.pi)

Dat levert na een kwartier wachten nog niets op. Ik bereken even de gravitatiepotentiaal omdat ik dan niet zit met de problematiek van het vectorieel karakter.

Als ik de absolute waarde in de noemer verander door een kwadraat, dan krijg ik waarschuwingen dat het niet goed convergeert.

EDIT: ik zie juist dat die absolute waarde een vierkantswortel moet zijn, maar voor het convergeren maakt dat niet veel uit.

[wnvl1]

Code: Selecteer alles

import numpy as np
from scipy.integrate import tplquad

a=6356000 #rstraal bij pool
b=6378000 #rstraal bij evenaar
c=6378000 #rstraal bij evenaar

G=6.6743e-11 

r_loc=1 #positie Chimay
phi_loc=0 #positie Chimay
theta_loc=np.radians(40) #positie Chimay

rho = 5515 #massa dichtheid aarde

tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc) )**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc)  )**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc)  )**2 )**0.5,
                0, 1,
                lambda theta: 0, np.pi,
                lambda theta, phi: 0, 2*np.pi)

Zo is beter.

[wnvl1]

Door de tolerantie te verlagen, komt er wel een uitkomst in een redelijke tijd.

Code: Selecteer alles

import numpy as np
from scipy.integrate import tplquad

a=6356000 #rstraal bij pool
b=6378000 #rstraal bij evenaar
c=6378000 #rstraal bij evenaar

G=6.6743e-11 

r_loc=1 #positie Chimay
phi_loc=0 #positie Chimay
theta_loc=np.radians(40) #positie Chimay

rho = 5515 #massa dichtheid aarde

tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc) )**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc)  )**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc)  )**2 )**0.5,
                0, 1,
                lambda theta: 0, np.pi,
                lambda theta, phi: 0, 2*np.pi,
       epsabs=0.0001, epsrel=0.0001)

(84552977.13936624, 26644.484128986347)

[wnvl1]

Voor de gravitatie kracht moet je dan een gelijkaardige integraal berekenen in de x, y en z richting (3 integralen dus) en dan de grootte berekenen van de 3 componenten tezamen. Of zou er een slimmere methode zijn?

[HansH]

als ik dat een maal heb dan kan mathcad dat heel snel uitrekenen

Of het sneller is dan mijn methode betwijfel ik.

qua methode zal het niet veel schelen, maar qua onderliggende engine voor de berekeningen wel. als die in C geschreven zijn kan het wel 1000 x zo snel zijn als een interpreter taal zoals python. maar snel is niet mijn eerste prioriteit hoewel het wel meegenomen is als het snel is.

[HansH]

ik ben er al uit:

er zit nog een foutje in de berekening, maar kan het bericht niet meer wijzigen, dus dan maar opnieuw posten.

er stonden 2 grenzen op -1 en 1 ipv op de ellipsoide gerelateerde waardes
berekende volume van de aarde klopt nu precies met de getallen hier: https://www.kuuke.nl/de-grootte-van-de-aarde/
dus moet dan wel goed zijn.

[wnvl1]

Als ik de gravitatiekracht projecteer op de x, y en z-as dan wordt het ondanks een verlaging van de tolerantie toch heel traag en loopt het geheugen vol (12GB). Het zijn natuurlijk ook complexe integralen.

Code: Selecteer alles

import numpy as np
from scipy.integrate import tplquad

a=6356000 #straal bij pool
b=6378000 #straal bij evenaar
c=6378000 #straal bij evenaar

G=6.6743e-11 

r_loc=1 #positie Chimay
phi_loc=0 #positie Chimay
theta_loc=np.radians(40) #positie Chimay

rho = 5515 #massa dichtheid aarde

#gravitatie: x-component
gravx = tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc) )**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc)  )**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc)  )**2 ) \
                *(r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))**2 + (r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc) )**2 + (r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc)  )**2 ) ,
                0, 1,
                lambda theta: 0, np.pi,
                lambda theta, phi: 0, 2*np.pi,
       epsabs=0.001, epsrel=0.001)

#gravitatie: y-component
gravy = tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc) )**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc)  )**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc)  )**2 ) \
                *(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc))/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))**2 + (r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc) )**2 + (r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc)  )**2 ) ,
                0, 1,
                lambda theta: 0, np.pi,
                lambda theta, phi: 0, 2*np.pi,
       epsabs=0.001, epsrel=0.001)

#gravitatie: z-component
gravz = tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc) )**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc)  )**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc)  )**2 ) \
                *(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc) )/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))**2 + (r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc) )**2 + (r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc)  )**2 ) ,
                0, 1,
                lambda theta: 0, np.pi,
                lambda theta, phi: 0, 2*np.pi,
       epsabs=0.001, epsrel=0.001)

[wnvl1]

De x en y component lukken bij heel hoge epsrel, op de z component loopt scipy vast.

Code: Selecteer alles

import numpy as np
from scipy.integrate import tplquad

a=6356000 #straal bij pool
b=6378000 #straal bij evenaar
c=6378000 #straal bij evenaar

G=6.6743e-11 

r_loc=1 #positie Chimay
phi_loc=0 #positie Chimay
theta_loc=np.radians(40) #positie Chimay

rho = 5515 #massa dichtheid aarde

#gravitatie: x-component
gravx = tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc))**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc))**2 )**1.5 *(r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))  ,
            0, 1,
            lambda theta: 0, np.pi,
            lambda theta, phi: 0, 2*np.pi,
            epsabs=0.1, epsrel=0.1)

print('Gravitatie x: ' +  str(gravx[0]) )

#gravitatie: y-component
gravy = tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc))**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc))**2 )**1.5 *(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc))  ,
            0, 1,
            lambda theta: 0, np.pi,
            lambda theta, phi: 0, 2*np.pi,
            epsabs=0.1, epsrel=0.1)

print('Gravitatie y: ' +  str(gravy[0]) )

#gravitatie: z-component
gravz = tplquad(lambda phi, theta, r: rho*r**2 *a*b*c* np.sin(theta)*G*1/((r*a*np.sin(phi)*np.cos(theta) -  r_loc*a*np.sin(phi_loc)*np.cos(theta_loc))**2 +(r*b*np.sin(phi)*np.sin(theta) - r_loc*b*np.sin(phi_loc)*np.sin(theta_loc))**2+(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc))**2 )**1.5 *(r*c*np.cos(phi) - r_loc*c*np.cos(phi_loc) )  ,
            0, 1,
            lambda theta: 0, np.pi,
            lambda theta, phi: 0, 2*np.pi,
            epsabs=0.1, epsrel=0.1)

print('Gravitatie z: ' +  str(gravz[0]) )

print((gravx[0]**2+gravy[0]**2+gravz[0]**2)**0.5)

[HansH]

met mathcad werkt het na wat hobbels opgelost te hebben tot nu toe goed.
ik kan het volume en de massa van de aarde berekenen op basis van integratie van volumedeeltjes met lokale dichtheid sigma.voor de sigma ben ik uitgegaan van de gegevens van xilvo met even simpel 1 rechte lijn gebruikt met waarde zodat de totale massa van de aarde klopt.

[Xilvo]

qua methode zal het niet veel schelen, maar qua onderliggende engine voor de berekeningen wel. als die in C geschreven zijn kan het wel 1000 x zo snel zijn als een interpreter taal zoals python. maar snel is niet mijn eerste prioriteit hoewel het wel meegenomen is als het snel is.

Ik gebruik voor de essentiële routine een JIT-compiler, dan is Python vrijwel net zo snel als gecompileerd C.
Dat is ook wel nodig, het scheelt ongeveer een factor 600 in snelheid.

Het grafiekje dat ik eerder plaatste (19 berekende punten met ruim 2 miljard massapunten per punt) kostte minder dan 9 minuten om te berekenen, zonder JIT-compiler was dat bijna 4 dagen geweest.

[Xilvo]

Voor de gravitatie kracht moet je dan een gelijkaardige integraal berekenen in de x, y en z richting (3 integralen dus) en dan de grootte berekenen van de 3 componenten tezamen. Of zou er een slimmere methode zijn?

Je zou meteen alleen de component in de richting van het gekozen punt aan het aardoppervlak kunnen berekenen.
Of dat slimmer, handiger en/of sneller is zou ik niet durven zeggen.

[wnvl1]

Kan je de richting waarin de gravitatie gaat staan op het aardoppervlak beredeneren vanuit symmetrie-overwegingen (uitgaande van een ellipsoïde)?