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?