Zaimplementowałem solver z Eulerem wstecznym w Pythonie 3 (używając numpy). Dla własnej wygody i jako ćwiczenie napisałem również małą funkcję, która oblicza przybliżoną różnicę skończoną różnicy gradientu, aby nie zawsze musiałem określać analitycznie jakobian (jeśli to w ogóle możliwe!).
Korzystając z opisów podanych w Ascher i Petzold 1998 , napisałem tę funkcję, która określa gradient w danym punkcie x:
def jacobian(f,x,d=4):
'''computes the gradient (Jacobian) at a point for a multivariate function.
f: function for which the gradient is to be computed
x: position vector of the point for which the gradient is to be computed
d: parameter to determine perturbation value eps, where eps = 10^(-d).
See Ascher und Petzold 1998 p.54'''
x = x.astype(np.float64,copy=False)
n = np.size(x)
t = 1 # Placeholder for the time step
jac = np.zeros([n,n])
eps = 10**(-d)
for j in np.arange(0,n):
yhat = x.copy()
ytilde = x.copy()
yhat[j] = yhat[j]+eps
ytilde[j] = ytilde[j]-eps
jac[:,j] = 1/(2*eps)*(f(t,yhat)-f(t,ytilde))
return jac
Przetestowałem tę funkcję, biorąc funkcję wielowymiarową dla wahadła i porównując symboliczny jakobian z gradientem wyznaczonym numerycznie dla zakresu punktów. Byłem zadowolony z wyników testu, błąd wynosił około 1e-10. Kiedy rozwiązałem ODE dla wahadła za pomocą przybliżonego jakobianu, zadziałało to również bardzo dobrze; Nie mogłem wykryć żadnej różnicy między nimi.
Następnie spróbowałem przetestować go za pomocą następującego PDE (równanie Fishera w 1D):
stosując dyskretyzację różnic skończonych.
Teraz metoda Newtona pojawia się w pierwszej chwili:
/home/sfbosch/Fisher-Equation.py:40: RuntimeWarning: overflow encountered in multiply
du = (k/(h**2))*np.dot(K,u) + lmbda*(u*(C-u))
./newton.py:31: RuntimeWarning: invalid value encountered in subtract
jac[:,j] = 1/(2*eps)*(f(t,yhut)-f(t,yschlange))
Traceback (most recent call last):
File "/home/sfbosch/Fisher-Equation.py", line 104, in <module>
fisher1d(ts,dt,h,L,k,C,lmbda)
File "/home/sfbosch/Fisher-Equation.py", line 64, in fisher1d
t,xl = euler.implizit(fisherode,ts,u0,dt)
File "./euler.py", line 47, in implizit
yi = nt.newton(g,y,maxiter,tol,Jg)
File "./newton.py", line 54, in newton
dx = la.solve(A,b)
File "/usr/lib64/python3.3/site-packages/scipy/linalg/basic.py", line 73, in solve
a1, b1 = map(np.asarray_chkfinite,(a,b))
File "/usr/lib64/python3.3/site-packages/numpy/lib/function_base.py", line 613, in asarray_chkfinite
"array must not contain infs or NaNs")
ValueError: array must not contain infs or NaNs
Dzieje się tak w przypadku różnych wartości eps, ale o dziwo tylko wtedy, gdy rozmiar kroku przestrzennego PDE i rozmiar kroku czasowego są ustawione tak, że warunek Courant-Friedrichs-Lewy nie jest spełniony. W przeciwnym razie to działa. (Jest to zachowanie, którego można się spodziewać, jeśli rozwiążesz go za pomocą Eulera!)
Dla kompletności, oto funkcja metody Newtona:
def newton(f,x0,maxiter=160,tol=1e-4,jac=jacobian):
'''Newton's Method.
f: function to be evaluated
x0: initial value for the iteration
maxiter: maximum number of iterations (default 160)
tol: error tolerance (default 1e-4)
jac: the gradient function (Jacobian) where jac(fun,x)'''
x = x0
err = tol + 1
k = 0
t = 1 # Placeholder for the time step
while err > tol and k < maxiter:
A = jac(f,x)
b = -f(t,x)
dx = la.solve(A,b)
x = x + dx
k = k + 1
err = np.linalg.norm(dx)
if k >= maxiter:
print("Maxiter reached. Result may be inaccurate.")
print("k = %d" % k)
return x
(Funkcja la.solve to scipy.linalg.solve.)
Jestem pewien, że moja poprzednia implementacja Eulera jest w porządku, ponieważ przetestowałem ją przy użyciu funkcji dla jakobianów i uzyskałem stabilne wyniki.
W debuggerze widzę, że newton () zarządza 35 iteracjami przed wystąpieniem błędu. Liczba ta pozostaje taka sama dla każdego eps, którego próbowałem.
Dodatkowa obserwacja: kiedy obliczam gradient za pomocą FDA i funkcji wykorzystując warunek początkowy jako dane wejściowe i porównuję oba, zmieniając jednocześnie wielkość epsilonu, błąd rośnie wraz ze zmniejszaniem się epsilonu. Spodziewałbym się, że początkowo będzie on duży, a następnie zmniejszy się, a potem znów będzie większy, gdy zmniejszy się epsilon. Zatem błąd w mojej implementacji jakobianów jest rozsądnym założeniem, ale jeśli tak, to jest tak subtelny, że nie jestem w stanie go dostrzec. EDYCJA: Zmodyfikowałem jacobian (), aby używać różnic naprzód zamiast różnic centralnych, a teraz obserwuję oczekiwany rozwój błędu. Jednak newton () nadal nie jest zbieżny. Obserwując dx w iteracji Newtona, widzę, że rośnie tylko, nie ma nawet wahań: prawie podwaja się (współczynnik 1,9) z każdym krokiem, przy czym współczynnik staje się coraz większy.
Ascher i Petzold wspominają, że przybliżenia różnic dla jakobianów nie zawsze działają dobrze. Czy przybliżony jakobian ze skończonymi różnicami może powodować niestabilność metody Newtona? A może przyczyną jest gdzie indziej? Jak inaczej mogę podejść do tego problemu?
Odpowiedzi:
Bardziej długi komentarz niż cokolwiek innego:
Spójrz na kod przybliżenia ilorazu różnicy SUNDIALS, aby uzyskać lepsze pojęcie o tym, co powinieneś zrobić w implementacji. Ascher i Petzold to dobra książka na początek, ale SUNDIALS jest faktycznie wykorzystywany w pracach produkcyjnych i dlatego został lepiej przetestowany. (Również SUNDIALS jest powiązany z DASPK, nad którym pracował Petzold.)
Empirycznie, przybliżone jakobianie mogą powodować niepowodzenia konwergencji w metodzie Newtona. Nie wiem, czy scharakteryzuję je jako „niestabilność”; w niektórych przypadkach po prostu nie jest możliwe osiągnięcie pożądanych tolerancji błędów w kryteriach zakończenia. W innych przypadkach może to objawiać się niestabilnością. Jestem prawie pewien, że zjawisko to jest bardziej ilościowe w książce metod numerycznych Highama lub w dyskusji Hairera i Wannera na temat metod W.
Zależy, gdzie według ciebie może być błąd. Jeśli jesteś bardzo pewny swojej implementacji wstecznego Eulera, nie zacznę od tego. Doświadczenie doprowadziło mnie do paranoi we wdrażaniu metod numerycznych, więc gdybym to był ja, zacznę od kodowania kilku naprawdę podstawowych problemów testowych (kilka niesztywnych i sztywnych problemów liniowych, równanie cieplne z wyśrodkowanym przybliżeniem skończonej różnicy, takie rzeczy) i użyłbym metody wytworzonych rozwiązań, aby upewnić się, że wiem, jakie będzie rozwiązanie i z czym powinienem się porównywać.
Jednak już to zrobiłeś:
To byłaby kolejna rzecz, którą przetestowałbym: użyj analitycznego jakobianu. Następnie możesz również przyjrzeć się ekstremalnym wartościom własnym skończonej różnicy jakobianów, gdy jesteś w niestabilnym regionie zacofanego Eulera. Spojrzenie na ekstremalne wartości własne analitycznego jakobianu jako podstawa do porównania może dać ci pewien wgląd. Zakładając, że wszyscy to sprawdzą, problem prawdopodobnie leży w rozwiązaniu Newtona.
źródło