18. Energia stanu podstawowego kationu cząsteczki wodoru z wykorzystaniem metody Monte Carlo i algorytmu Metropolisa

Przygotowany uprzednio program do obliczania energii stanu podstawowego atomu wodoru z wykorzystaniem metody Monte Carlo i algorytmu Metropolisa zaadaptowałem do obliczenia energii stanu podstawowego cząsteczki H₂⁺. Funkcja falowa, którą wybrałem do obliczeń była kombinacją liniową orbitali 1s atomów wodoru 1:

Parametrem obliczeniowym - odległość między protonami R_j w układzie. Rysunek poniżej przedstawia sytuację geometryczną użytą do obliczeń:

W wyniku przeprowadzonych obliczeń otrzymałem zbiór danych, które są przedstawione na poniższym wykresie:

Wartość minimalna energii wynosi około -0.567 j.at. przy odległości między protonami wynoszącej 2.5 angstrema.
Tekst programu, który wygenerował zbiór wyników jest przytoczony poniżej:

DECLARE FUNCTION fufal! (x!, y!, z!)
DECLARE FUNCTION GRNF! ()  'gaussowski generator pseudolosowy
                           'jednostki atomowe
llos = 500000            'liczba losowan
h = .01                  'krok rozniczkowania
h2 = h * h               'kwadrat kroku rozniczkowania
DIM SHARED Rj            'odl. miedzy protonami
OPEN "H2plus.txt" FOR APPEND AS #1
FOR Rj = .9 TO 6 STEP .1    'obl. energi w zal. od Rj

En = 0                   'zerowanie energii do sumowania
x = 0: y = 1.2: z = 0    'poczatkowy punkt lancucha Markowa
FOR i = 1 TO llos
  xzapas = x             'zabezpieczenie starych wspolrzednych
  yzapas = y
  zzapas = z

  psi = fufal(x, y, z)    'wartosc psi w punkcie (x,y,z)
  psikw = psi * psi       'kwadrat psi

           'Losowanie nowego punktu  (xn,yn,zn)
  xn = 2 * RND - 1        'zakres (-1, 1)
  yn = 2 * RND - 1
  zn = 2 * RND - 1

  pr = SQR(xn*xn + yn*yn + zn*zn)  'promien punktu (xn,yn,zn)
  skok = .6 * GRNF         'skok przypadkowy wg rozkladu Gaussa

  x = x + xn * skok / pr  'poprawione wspolrzedne nowego punktu
  y = y + yn * skok / pr
  z = z + zn * skok / pr
 
  psinowa = fufal(x, y, z)     'wartsc psi w nowym punkcie
  psinowakw = psinowa * psinowa

'obliczanie warunku przejscia lancucha do nowego punktu

 rel = psinowakw / psikw   'stosunek nowej i starej psi
  IF rel < 1 THEN          'przejscie bezwarunkowe gdy rel >= 1
    IF RND > rel THEN      'przejscie odrzucone
      x = xzapas
      y = yzapas
      z = zzapas
      psinowa = fufal(x, y, z)
    END IF
 END IF

 'obliczanie calki z energii lokalnych
  'obliczanie pochodnej w punkcie (x,y,z)
  fufal2 = 2 * fufal(x, y, z)
  d1 = (fufal(x - h, y, z) - fufal2 + fufal(x + h, y, z)) / h2
  d2 = (fufal(x, y - h, z) - fufal2 + fufal(x, y + h, z)) / h2
  d3 = (fufal(x, y, z - h) - fufal2 + fufal(x, y, z + h)) / h2
 
  Lappsi = -(d1 + d2 + d3) / 2'druga poch. psi w punkcie (x,y,z)

  V1 = -1/(SQR((x - Rj/2)^2 + y*y + z*z))'potencjal  psi/psi=1
  V2 = -1/(SQR((x + Rj/2)^2 + y*y + z*z))
  En = En + Lappsi/psinowa+V1+V2+1/Rj    'suma energii lok.
 NEXT i

   PRINT #1, Rj, En / llos   'druk zm. wariacyjnej i en. calk.
NEXT Rj
  CLOSE #1
END

FUNCTION fufal (x, y, z)
  x1 = x - Rj / 2
  x2 = x + Rj / 2
  r1 = SQR(x1 * x1 + y * y + z * z)
  r2 = SQR(x2 * x2 + y * y + z * z)
  fufal = EXP(-r1) + EXP(-r2)
END FUNCTION

FUNCTION GRNF
  r1 = SQR(-LOG(1 - RND))
  GRNF = r1 * SIN(2 * 3.14159265# * RND)
END FUNCTION

17. Dystrybucja zmiennej x dla kwadratu funkcji logarytmicznej w obszarze zachowania chaotycznego

W poprzednim poście pokazałem, że kwadrat funkcji logarytmicznej posiada wykładnik Lapunowa o wartościach powyżej zera. Poniżej przedstawiam wynik działania programu, który liczy dystrybucję zmiennej x w obszarze chaotycznym:

'Dystrybucja x dla kwadratu funkcji logarytmicznej
CLS : SCREEN 12
vec = 510    'rozmiar wektora zbierana danych
DIM d(vec)
x = .7
llos = 10000  'liczba iteracji
c = 1.2       'dla c=1.2 min=1.15E-9, max=507.92
OPEN "dystryb.txt" FOR APPEND AS #1
FOR i = 1 TO llos
  x = c * LOG(x) ^ 2  'funkcja logarytmiczna
  yd = INT(x)         'wskaznik w wektorze
  d(yd) = d(yd) + 1
NEXT i

FOR i = 1 TO vec
   PRINT #1, i, d(i)      'drukowanie wynikow
  'PSET (i, 400 - d(i))
NEXT i
CLOSE #1

Wynikiem działania programu jest zbiór danych, które przedstawiłem na poniższym wykresie:

Bez szczegółowej analizy można powiedzieć, że dystrybucja przypomina przebiegiem funkcję 1/x.

16. Wykładnik Lapunowa dla kwadratu funkcji logarytmicznej

Zaobserwowałem, że funkcja logarytmiczna może zachowywać się chaotycznie. W celu potwierdzenia takiej możliwości obliczyłem wykładniki Lapunowa dla szeregu wartości stałej c w równaniu 1:

Zadanie to realizuje program przedstawiony poniżej:

'wykladnik Lapunowa
'dla rownania logarytmicznego
CLS : SCREEN 12
N = 4000      ' liczba iteracji
OPEN "lapunow3.txt" FOR APPEND AS #1
FOR c = 1.1 TO 3 STEP .01    'wspolcz. w rown. logarytmicznym
  suma = 0
  x = .7
  FOR i = 1 TO N
    suma = suma + LOG(ABS(2 * c * LOG(x) / x))'suma log pochod.
    x = c * LOG(x) ^ 2 'iterowane rown. logarytmiczne
  NEXT i
  Lapunow = suma / N    ' wstpoczynnik Lapunowa
  PRINT #1, c, Lapunow
  'PSET (c * 150, 100 - Lapunow * 100)  'obraz na ekranie
NEXT c
CLOSE #1

Program generuje zbiór wartości, które są przedstawione na poniższym wykresie. Należy dodać, że wykres został skompilowany z dwóch zbiorów dla c (0.1,1.0) oraz (1.1, 3.0). Podzial był związany z błedami zaokrągleń zmiennej c. Zaokrąglenia powodowały, że pojawiły się błędy dla dla wartości c około 1.58.

Wykres pokazuje, że chaotyczne zachowanie kwadratu logarytmu jest możliwe dla c w obszarze od około 0.10 do około 1.84.

15. Dystrybucja zmiennej równania logistycznego dla obszaru chaotycznego

Przedmiotem tego wpisu jest chaotyczne zachowanie się równania logistycznego dla pewnych wartości stałej c danego wzorem 1:

Zadałem sobie pytanie, ile razy wybierana jest konkretna wartość x_i podczas procesu iteracji równania 1. Zagadnienie to rozwiązałem przy użyciu programu, który zlicza pojawiające się wartości zmiennej x a wynik zliczania wpisuje do wektora d. W celu "rozciągnięcia" zakresu zmiennych każda wartość x_i jest mnożona przez 10000. Tekst programu przedstawiłem poniżej:

'dystrybucja zm. x funkcji logistycznej
DIM d(10000)
y = .7
llos = 1000000
c = 3.8
OPEN "dystryb.txt" FOR APPEND AS #1
FOR i = 1 TO llos
  x = y
  y = c * x * (1 - x)
  yd = INT(y * 10000)
  d(yd) = d(yd) + 1
NEXT i

FOR i = 1 TO 10000
  PRINT #1, i, d(i)
NEXT i
CLOSE #1

Wynikiem działania programu jest zbiór danych, ktore zamieszczone są poniżej:

Można zauważyć, że średnio rzecz biorąc, najliczniej prezentowana wartość około 250. Liczności wzrastają na brzegach przedziału. Wspomnę, że zastosowałem pojedynczą precyzję w deklaracjach zmiennych. Zasadniczo powinienem użyć definicji zmiennych w podwójnej precyzji. Ten zabieg powinien spowodować "wygładzenie" poszarpanego wykresu dystrybucji.

14. Wykładnik Lapunowa dla rownania logistycznego. Wprawka I

Wykładnik Lapunowa jest dla mnie interesujący z tego względu, że pozwala określić wartości współczynników w równaniach, które mogą zachowywać się w sposób chaotyczny. Chcąc stwierdzić, czy dobrze rozumiem istotę zagadnienia napisałem program w języku BASIC, który liczy wykładniki Lapunowa dla równania logistycznego postaci 1:

Wykładnik Lapunowa określony jest wzorem 2:

przy czym każdą następną wartość x oblicza się ze wzoru 1.
Tekst programu znajduje się poniżej. Program umożliwia zapisanie danych na dysku lub obejrzenie wyniku w postaci wykresu bezpośrednio na ekranie:

'wykladnik Lapunowa
'dla rownania logistycznego
SCREEN 12
N = 4000      ' liczba iteracji
x = .7        ' punkt poczatkowy
CLS : SCREEN 12
OPEN "lapunow.txt" FOR APPEND AS #1
FOR c = 2.1 TO 3.9 STEP .001  'wspolcz. w rown. logistycznym
  suma = 0
  FOR i = 1 TO N
    suma = suma + LOG(ABS(c - 2 * c * x)) 'suma log pochodnych
    x = c * x * (1 - x)  'iterowane rown. logistyczne
  NEXT i
  Lapunow = suma / N     'wykladnik Lapunowa
  PRINT #1, c, Lapunow
  'PSET (c * 150 - 200, 200 - Lapunow * 100) 'obraz na ekranie
NEXT c
CLOSE

Wynik działania programu znajduje się poniżej w postaci wykresu:

Dla wartości wykładnika Lapunowa >0 równanie zachowuje się chaotycznie.