Symulacja rozłożenia ładunków na sferze

1

Kładziemy dwa jednakowe ładunki na sferze, wtedy one ustawią się na pewno po przeciwnych stronach - z powodu odpychania.

Bierzemy więcej ładunków, np. 6 sztuk - i jak one się teraz rozłożą?
Pewnie symetrycznie, po dwa na każdej z osi: x,y, z, po przeciwnych stronach: +/-r, czyli utworzą ośmiościan.

Ale co dalej, np. 9, albo 13,... 37 - nie ma przecież takiego wielościanu foremnego, na którym mogłyby się poustawiać symetrycznie?!

Więcej np. n = 100 - no i jak to się poustawia, aby zachować perfekcyjną symetrię na sferze?

A może to się nigdy nie ustabilizuje,
znaczy ładunki ustawione na sferze będą biegać po jej powierzchni w nieskończoność, szukając równowagi, która jest niemożliwa? :)

0

Wydaje mi się (z fizyką miałem ostatnio do czynienia na studiach), że po prostu wpadną w którąś ze studni potencjału — np. ustawią się na jednej płaszczyźnie na wielokącie foremnym.

0
Althorion napisał(a):

Wydaje mi się (z fizyką miałem ostatnio do czynienia na studiach), że po prostu wpadną w którąś ze studni potencjału — np. ustawią się na jednej płaszczyźnie na wielokącie foremnym.

Takie coś byłoby skrajnie nietrwałe.

Weź np. 32 ładunki i ustaw w kółeczku na tej sferze;
potem pyknij jeden lekko w bok, a on od razu poleci daleko i jak rakieta, bo z powodu silnego odpychania tego całego koła 31 ładunków!

0

Jak rozważasz ładunki na sferze w kompletnej izolacji od otoczenia, to nie będzie miało tego co „pyknąć”. A jak nie, to regularnie owo otoczenie będzie zaburzać układ tak czy owak.

0

Ok, Zatem taka symulacja powinna wyglądać następująco:

  1. ustawiamy na starcie 1 ładunek, więc on w dowolnym punkcie sfery pozostanie na miejscu
  2. dokładamy drugi, kolejny - w dowolnym miejscu sfery, więc ten już nie ustoi, bo będzie odpychany przez ten pierwszy;
    zacznie biegać po tej sferze, a ten pierwszy również - i symetrycznie, i nigdy się nie zatrzymają!
  3. zatem aby to zatrzymać wprowadzamy tarcie, aby mogły w końcu wyhamować,
    zależne od prędkości, np. tarcie: T = -kv^2, albo dowolna inna funkcja: T = -k exp(|v|), itp.

W sumie można od razu wstawić n ładunków w losowych pozycjach na sferze;
one zaczną wtedy tańczyć, ale to tarcie musi je w końcu zatrzymać;
no więc to będzie ta szukana konfiguracja, ale niekoniecznie jedyna możliwa.

0

Nie będą się po prostu układać jak pierwiastki n-tego stopnia z 1 na dowolnym okręgu bedacym przekrojem wspomnianej sfery? (Przekroj plaszczyzna przechodząca przez srodek)

1

W mechanice klasycznej taki układ musi mieć minimum energii, więc musi istnieć układ stabilny.
Jedyna równica jak może być, to że istnieje wiele nietrywialnie różnych minimów o tej samej energii potencjalnej.

Dodatkowo, jako że cały układ ma symetrię sferyczną, to każde minimum podlega obrotom, które nic nie zmieniają, więc nie ma minimum w postaci punktu w przestrzeni fazowej.

0
yarel napisał(a):

Nie będą się po prostu układać jak pierwiastki n-tego stopnia z 1 na dowolnym okręgu bedacym przekrojem wspomnianej sfery? (Przekroj plaszczyzna przechodząca przez srodek)

Tak byłoby w wersji 2D, czyli na okręgu, a tu masz sferę - wersja 3D.

2 i 3 ładunki będą na okręgu, bo tu nie ma wyjścia: każde 3 punkty leżą w płaszczyźnie!
No, a 4 już utworzą czworościan.

5 sztuk - no i jak to się ustawi?
Na pewno tak: 3 ustawią się w trójkąt, a 2 pozostałe będą skrajnie daleko, czyli na osi pionowej, prostopadłej do tego trójkąta.

6 - ośmiościan
7 - ?
8 - ? (nie będzie to sześcian, wbrew pozorom)
9 - ? pewnie 3 trójkąty.
...
12 - pewnie 20-ścian foremny.
14 - chyba jest wielościan
...

24 - tu jest chyba bryła symetryczna - do tego taka dualistyczna, czyli w dwóch wersjach...
60 - tu jest podobnie jak z 24

1
yarel napisał(a):

Nie będą się po prostu układać jak pierwiastki n-tego stopnia z 1 na dowolnym okręgu bedacym przekrojem wspomnianej sfery? (Przekroj plaszczyzna przechodząca przez srodek)

Masz na myśli pierwiastki zespolone tak?

Raczej nie rozłożą się na jednym okręgu, wciąż miałyby energię potencjalną wyższą niż w przypadku rozłożenia na całej kuli ;) Energia potencjalna ładunku maleje z kwadratem odległości, więc ogólnie będą dążyły do wyrównania odległości między każdą parą ładunków. Podejrzewam że dla uproszczenia można by wziąć jakąś projekcję i sprowadzić problem do zagadnienia 2D.

Idąc tym tokiem zdefiniowałem sobie funkcję celu (opartą o iloczyn skalarny - w końcu im bardziej przeciwne punkty, tym mniejszą da wartość, nie potrzebuję przecież liczyć energii potencjalnej ;) ) i minimalizujemy sumę iloczynów skalarnych (tu maksymalizujemy liczbę przeciwną i jeszcze przeliczoną żeby maksimum było równe 1, ale to szczegół):

import math
import numpy as np

def to_3D(array: np.ndarray):
    arr = []
    for pt in array[:]:
        arr.append([
            math.cos(pt[0]) * math.cos(pt[1]),
            math.sin(pt[0]) * math.cos(pt[1]),
            math.sin(pt[1])
        ])
    return np.array(arr)
        
def dist_3D_fun(array: np.ndarray):
    _3d_vecs = to_3D(array)
    result = _3d_vecs @ _3d_vecs.transpose()
    return 1.0 - np.sum(result)/result.shape[0]

Zakładam, że proponowane rozwiązania to kąty względem osi OX i OZ - bo tak ;) więc potrzebuję przeliczenia na współrzędne kartezjańskie, by dostać iloczyn skalarny - dla uproszczenia sfera ma zerowy promień.

No to teraz odpalamy sobie optymalizację, kiedyś tam z nudów naskrobałem na kolanie klasę Specimen z bieda-operatorami genetycznymi na np.ndarray'sach, więc je wykorzystam (tutaj dla dwóch punktów, stąd shape [2, 2]:

threshold = 1e-7
steps_limit = 10000
dimensions = [2,2]
a = Specimen(dimensions)
b = Specimen(dimensions)
a(dist_3D_fun)
b(dist_3D_fun)
c = a + b
steps = 0
delta = threshold
while delta >= threshold and steps < steps_limit:
    steps += 1
    if a.performance < b.performance:
        d = b
        b = a
        a = d
    # print("Step {0}: {1:.9f}".format(steps, a.performance))
    if c.performance >= a.performance:
        delta = c.performance - a.performance
        b = a
        a = c
    elif c.performance >= b.performance:
        delta = c.performance - b.performance
        b = c
    c = a + b
    c.mutate(1.0/(steps ** 0.5))(dist_3D_fun)
    
print("Step {0}: {1:.3f}".format(steps, a.performance))
print(to_3D(a._genes))

No i wskazało przeciwległe punkty (z jakimś tam niewielkim błędem na 3-cim/4-tym miejscu po przecinku):

Step 4228: 1.000
[[-0.59341103  0.4789695  -0.64687832]
 [ 0.59333623 -0.47878774  0.64708147]]

No to testujemy dla jakiejś ciekawszej liczby, np. n=11 (shape [11, 2]):

Step 10000: 0.060
[[ 0.92665822  0.36947295  0.06924074]
 [-0.15082731 -0.28915089  0.94532687]
 [ 0.94456119 -0.2580341   0.20303339]
 [-0.5927377  -0.0256576   0.80498678]
 [ 0.57147413  0.18094127  0.80042337]
 [ 0.73409681  0.34770983 -0.58326645]
 [-0.81766341  0.03907948  0.57436865]
 [-0.07086159 -0.23513502  0.96937617]
 [-0.33180022 -0.1300778   0.93433847]
 [-0.00689443  0.00534034 -0.99996197]
 [ 0.36834138  0.14222481 -0.91874737]]

Wynik 0.060 brzmi słabo, jak tak się przyjrzymy to część kropek faktycznie jest skupiona w jednym miejscu (choć nie ogarnąłem jeszcze, jak wyświetlić interaktywny wykres 3D w Jupyter Notebook):
screenshot-20181118234702.png

Próbujemy jeszcze raz:

Step 5990: 1.000
[[ 0.23065809 -0.07382626  0.97023014]
 [-0.29096953  0.77687333 -0.55839463]
 [ 0.52335038 -0.81229086 -0.25746445]
 [-0.49236392 -0.30186805 -0.816366  ]
 [ 0.13249796 -0.65986807 -0.73960694]
 [-0.14103619  0.03070679 -0.98952811]
 [-0.39512705  0.24833482  0.88442322]
 [ 0.11671746  0.22656547 -0.96697731]
 [-0.15520067  0.43405139  0.88741881]
 [ 0.76863071  0.10062412  0.63172907]
 [-0.29581929  0.02816136  0.95482872]]

I dostajemy wykres - na pierwszy rzut oka wygląda to trochę lepiej, choć ciężko stwierdzić nie mogąc sobie poobracać i poprzeglądać z każdej strony jak te ładunki nam się rozłożyły:
screenshot-20181118234928.png

@exp co do ustawienia pięciu, dostałem przykładowo coś takiego:

Step 3749: 1.000
[[-0.47633634 -0.04123628  0.87829566]
 [-0.2678149   0.39115383 -0.88049637]
 [-0.85745299  0.50291791 -0.10884828]
 [ 0.67367885 -0.70312072 -0.22754794]
 [ 0.92800035 -0.15393657  0.33929174]]

screenshot-20181118235312.png

0

Właściwie to skoro dla N punktów na N-1 wymiarowej sferze rozwiązaniem zawsze będzie simpleks wpisany w tę sferę (trójkąt równoboczny w okręgu, czworościan foremny w sferze itd) to wystarczyłoby w sumie obliczyć współrzędne tych punktów na N-1 wymiarowej sferze i potem zredukować liczbę wymiarów do 3, ale nie mam pomysłu jak to zrobić, obcinanie nadmiarowych współrzędnych zwróci pewnie głupoty.

0

Iloczyn skalany raczej nie za bardzo nadaje się do tego.

Cyba tylko dla dwóch wychodzi dobrze, bo: cos180 = -1 jest tu minimum;

A co dla trzech wyjdzie?
Chyba też -1, czyli trójkąt prostokątny: -1 + 0 + 0 = -1

Dla czterech pewnie kwadrat: -1 - 1 + 0 + 0 + 0 ... = -2
bo dla czworościanu byłoby to większe - chyba dodatnie.

Zatem lepiej użyć odpychania w stylu sprężyny, co daje potencjał: -kr^2,
więc zamiast iloczynu skalarnego położeń: r1*r2, wystarczy tu wyliczyć wpierw różnicę: r = r1-r2 = (x1-x2, y1-y2, z1-z2),
i wymnożyć to teraz ale ze sobą i zanegować wynik: -r.r

co niewiele zmienia algorytm, a wyniki będą znacznie bardziej sensowne.

0
exp napisał(a):

Iloczyn skalany raczej nie za bardzo nadaje się do tego.

Cyba tylko dla dwóch wychodzi dobrze, bo: cos180 = -1 jest tu minimum;

A co dla trzech wyjdzie?
Chyba też -1, czyli trójkąt prostokątny: -1 + 0 + 0 = -1

Raczej równoboczny: -1/2 + -1/2 + -1/2 = -3/2

Dla czterech pewnie kwadrat: -1 - 1 + 0 + 0 + 0 ... = -2
bo dla czworościanu byłoby to większe - chyba dodatnie.

Dla czworościanu również wychodzi -2 (6 * -1/3). Nie może być dodatnie, kąty nadal będą większe od PI/2.

Zatem lepiej użyć odpychania w stylu sprężyny, co daje potencjał: -kr^2,
więc zamiast iloczynu skalarnego położeń: r1*r2, wystarczy tu wyliczyć wpierw różnicę: r = r1-r2 = (x1-x2, y1-y2, z1-z2),
i wymnożyć to teraz ale ze sobą i zanegować wynik: -r.r

co niewiele zmienia algorytm, a wyniki będą znacznie bardziej sensowne.

Wydaje mi się, że różnicy to wielkiej nie zrobi, i tak funkcja będzie mogła wpaść w minimum lokalne.

0

Tak, miałem na myśli pierwiastki zespolone z jedynki. Czy odległość na sferze liczona jest z uwzględnieniem krzywizny sfery, czy po prostu jako odległość euklidesowa w 3D ?

Ja bym jeszcze spojrzał przez współrzędne sferyczne, tj. punkt na sferze identyfikujemy przez trójkę (R,Theta, Phi) - promień sfery i 2 kąty.

Dla N punktów:

  1. Losujemy kąty Theta, Phi z przedziału [0,2pi]
  2. Generujemy Theta + 2pi*k/N, dla k=0..N-1
  3. Generujemy Phi+2pi*j/N, dla j=0..N-1
  4. Robimy kartezjana i dostajemy N^2 punktów (współrzędne kątowe, R olewamy bo nie wnosi nic ciekawego do tematu)
  5. Uruchamiamy K-means na takim N^2 zbiorze punktów i szukamy N centroidów

Nie wiem czy taki porządny rozkład punktów w 4) nie wykrzaczy K-means ;-)

0

Szukałem wątku o nauce pod Data Science, ale nie wytrzymałem i się zarejestrowałem ;).

Moje 3 grosze:
rozkład musi minimalizować potencjał. Dla dyskretnych ładunków potencjał jest dany równaniem (pisza tak na wpol LaTex-em)
Potencjał P P=  \sum_{i,j, i \neq j} \frac{1}{|\vec r_i - \vec r_j|}
gdzie \vec r_i = \cos(a_i) \sin(b_i) \hat x + \sin(a_i) \sin(b_i) \hat y + \cos(b_i) \hat z
gdzie \hat x to wektor w kierunku x.
% sumujemy po indeksie ładunku. W sumie unikamy sytuacji i == j , tzw. samooddziaływania. Powoduje ono ekspolozję kosmosu tzn. dzielimy przez zero i otrzymujemy nieskonczoną energię.*

Trzeba znaleźć minimum takiego potencjału. Żeby sobie sprawę ułatwić można założyć że jeden z punktów zawsze bedzie na szczycie kuli \vec r_1 = 1 \hat z (zakładając że kula ma promien r=1). To załatwi część problemu jaki stanowi symetria - każdy rozkład obrócony o kąt wokół dowolnej osi da równie dobre rozwiązanie. Jak "przymocujemy" jeden punkt na osi Z to zostaje nam tylko problem obrotu wokół osi Z. Dlatego drugi punkt też mocujemy, ale tylko poprzez ustalenie na jakim kole leży. I tak np \vec r_2 =  sin(b_i) \hat x  + cos(b_i) \hat z (leży na kole rozpiętym na osiach X i Z).

Mam nadzieję, że mój wywód jest zrozumiały. Fizykę umiem (robię PhD z tej pięknej dziedziny), ale nie znam się dobrze na metodach numerycznych (trzeba znaleźć minimum tej funkcji, o ile się nie mylę fizyka gwarantuje, że powinno istnieć tylko jedno minimum będące tym samym minimum globalnym).

Jeżeli coś jest niezrozumiałe to chętnie rozwinę.

*W prostej elektrodynamice macha się rękami i mówie, że to jest powód dla którego istnieje prawo zachowania ładunku. Energia związana z samoodziaływaniem kojarzona jest z energią potrzebną do stworzenia ładunku. Ładunku stworzyć nie można (bo potrzeba nieskończenie wiele energii) więc ładunek zostaje choćby nie wiem co. Nawet czarne dziury muszą się podporządkować.

0
superdurszlak napisał(a):
exp napisał(a):

Iloczyn skalany raczej nie za bardzo nadaje się do tego.

Cyba tylko dla dwóch wychodzi dobrze, bo: cos180 = -1 jest tu minimum;

A co dla trzech wyjdzie?
Chyba też -1, czyli trójkąt prostokątny: -1 + 0 + 0 = -1

Raczej równoboczny: -1/2 + -1/2 + -1/2 = -3/2

zgada się: 3cos120 = -3/2

Dla czterech pewnie kwadrat: -1 - 1 + 0 + 0 + 0 ... = -2
bo dla czworościanu byłoby to większe - chyba dodatnie.

Dla czworościanu również wychodzi -2 (6 * -1/3). Nie może być dodatnie, kąty nadal będą większe od PI/2.

w czworościanie tak będzie: (1,1,1)(-1,-1,1)/3 = -1/3,
i jest 6 takich kombinacji:
6* -1/3 = -2, czyli tyle samo co dla kwadratu.

Wydaje mi się, że różnicy to wielkiej nie zrobi, i tak funkcja będzie mogła wpaść w minimum lokalne.

W przypadku potencjału 1/r, albo r^2 wyniki byłby raczej inne - można sprawdzić ten czworościan...

sprawdźmy co to siły mają potencjał w stylu: V = r1.r2 = cos(f);

co chyba daje siłę typu:
F = sin(f)/r * f^0 = r1 x r2 / r, iloczyn wektorowy...

0

Sprawdźmy co tu wyjdzie dla n = 4 w tych dwóch wersjach czworościan i kwadrat;

A. z potencjałem 1/r (bez minusa bo tu jest odpychanie):

odległości są jednakowe - bok czworościanu: a^2 = 8/3 (dla R=1)
zatem potencjał: V = 6/a = 6 sqrt(3/8) = 3.674

a teraz kwadrat: V = 1/2 + 2*1/sqrt2 + 1/2 + 1/sqrt2 + 1/sqrt2 = 1 + 4/sqrt2 = 3.828

zatem czworościan ma mniejszy potencjał, więc wygrywa.

B. wersja z iloczynem skalarnym: V = ri.rj;

tu otrzymamy -2 i -2, czyli jednakowo.

C. jeszcze potencjał sprężynowy, czyli typu: r^2 (tradycyjne: F = kx);

czworościan: 6 * a^2 = 6 * 8/3 = 16
kwadrat: 4 + 2*2 + 4 + 2 + 2 = 16

czyli to samo... bardzo dziwne - przecież kwadrat nie jest stabilny, więc jakim cudem może być taki sam potencjał?

Prawdopodobnie w wersji z potencjałem: V = r^2 istnieje lepsze rozwiązanie - pośrednie pomiędzy kwadratem i czworościanem.

0
AndrzejGajos napisał(a):

*W prostej elektrodynamice macha się rękami i mówie, że to jest powód dla którego istnieje prawo zachowania ładunku. Energia związana z samoodziaływaniem kojarzona jest z energią potrzebną do stworzenia ładunku. Ładunku stworzyć nie można (bo potrzeba nieskończenie wiele energii) więc ładunek zostaje choćby nie wiem co. Nawet czarne dziury muszą się podporządkować.

Energia ładunku Q rozłożonego równomiernie na sferze wynosi: 1/2 Q^2/R, zatem nie ma z tym raczej problemów - R > 0.

Tyle że takie coś wymaga ekstra sił z zewnątrz do utrzymania - ciśnienia, no i stąd pomysły z wirującymi elektronami - wiry, spiny, magnetyzm ściskający, itp. bajki.

1 użytkowników online, w tym zalogowanych: 0, gości: 1