F1(N) - stablicowana funkcja (dokładniej wygenerowany, znormalizowany rozkład gaussa N(0,4)), N=500
F2(N) - funkcja F2=F1'+constF1=\(\displaystyle{ \frac{F(i-1)-2F(i)+F(i+1)}{h^2}+CF1}\) (trójpunktowy wzór na drugą pochodną)
Teraz zapiszę skrótowo kod, z którym jest problem w "pseudo-fortranie":
Kod: Zaznacz cały
PROGRAM
DO i=1,N !pętla po wszystkich elementach (wartościach funkcji)
DO t=0.0,Kdt,dt !pętla "czasowa" - krok dt, K iteracji
F1=F1-F2 !olbrzymi skrót - tu oczywiście korzystam z procedury, by dokonać czasowej ewolucji funkcji F1, w tym celu również ewoluuję F2
END DO
WRITE F1 !wypisanie "wyewoluowanej" F1 po K iteracjach
END DOMoje pytania to oczywiście: dlaczego i jak sobie z tym poradzić?
Moje badanie:
F2 z powodu małego kroku h jak i małego licznika mogłoby powodować ≈ 0/0, a więc symbol nieoznaczony, czyli NaN. Ten problem też się pojawił, ale go wyeliminowałem - teraz F2 - przed wywołaniem procedury ewolucji sprawuje się wyśmienicie (więc warunki początkowe są OK).
SUBROUTINE: F1=F1-F2 - dla pewnych wartości działa dobrze (np. gdy ewoluuję F1 w x=0), nie powinno też generować żadnych NaN-ów, bo to tylko różnica dwóch małych liczb, a więc najwyżej powinienem dostać ładne zaokrąglone 0.
Problem siedzi więc (chyba) w dalszych ewolucjach F2, dla naprawdę małych wartości F1 (trochę dalej od maksimum krzywej w x=0)
Jakby to ciachnąć?
Problemem jest zapewne ograniczona precyzja liczb zmiennoprzecinkowych - może zwiększyć REAL na DOUBLE PRECISION? (choć sprawdzałem wcześniej jak to się zmienia w różnych przypadkach i mała poprawa...)
Może zmniejszyć przedział obliczeń? Prowadzę je dla x=[-250;250] - wiem, że to straszne biorąc pod uwagę funkcję jaką się zajmuję, ale takie były jednak odgórne polecenia (btw, h=1.0).
Mój faworyt: może przeskalować tego Gaussa - wartości x 100 np.?
Help.
Zaznaczam, powyżej to pseudokod (nieprecyzyjny - jedynie obrazowy) - wszystko się dobrze kompiluje, nie ma warnów, szczególiki są sprawdzone, warunki brzegowe (zerowanie F1 i F2 dla granic przedziału) i początkowe dobre, i poprawnie nałożone.
Jeżeli ktoś się zainteresuje i będzie chciał zobaczyć jak to naprawdę wygląda, to oczywiście udostępnię kod. Nie zrobiłem tego teraz, żeby nie straszyć - trochę tego jest, trochę skomplikowane, a ponadto z pewnością dalekie od programistycznej estetyki.
