J'ai implémenté un solveur Euler vers l'arrière en python 3 (en utilisant numpy). Pour ma commodité et comme exercice, j'ai également écrit une petite fonction qui calcule une approximation par différence finie du gradient afin que je n'ai pas toujours à déterminer le jacobien analytiquement (si c'est même possible!).
En utilisant les descriptions fournies par Ascher et Petzold 1998 , j'ai écrit cette fonction qui détermine le gradient à un point donné 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
J'ai testé cette fonction en prenant une fonction multivariée pour le pendule et en comparant le jacobien symbolique au gradient déterminé numériquement pour une gamme de points. J'ai été satisfait des résultats du test, l'erreur était d'environ 1e-10. Lorsque j'ai résolu l'ODE pour le pendule en utilisant le jacobien approximatif, cela a également très bien fonctionné; Je n'ai pu détecter aucune différence entre les deux.
Ensuite, j'ai essayé de le tester avec la PDE suivante (équation de Fisher en 1D):
en utilisant une discrétisation par différence finie.
Maintenant, la méthode de Newton explose au premier pas de temps:
/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
Cela se produit pour une variété de valeurs eps, mais étrangement, uniquement lorsque la taille du pas spatial PDE et la taille du pas de temps sont définies de sorte que la condition Courant – Friedrichs – Lewy n'est pas remplie. Sinon ça marche. (C'est le comportement que vous attendez si vous résolvez avec Euler avant!)
Pour être complet, voici la fonction de la méthode Newton:
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
(La fonction la.solve est scipy.linalg.solve.)
Je suis convaincu que mon implémentation en arrière d'Euler est en ordre, car je l'ai testée en utilisant une fonction pour le jacobien et obtenir des résultats stables.
Je peux voir dans le débogueur que newton () gère 35 itérations avant que l'erreur ne se produise. Ce nombre reste le même pour chaque eps que j'ai essayé.
Une observation supplémentaire: lorsque je calcule le gradient avec FDA et une fonction en utilisant la condition initiale comme entrée et que je compare les deux tout en faisant varier la taille d'Epsilon, l'erreur augmente à mesure qu'Epsilon se rétrécit. Je m'attendrais à ce qu'il soit grand au début, puis devienne plus petit, puis à nouveau plus grand à mesure que epsilon se rétrécit. Une erreur dans ma mise en œuvre du jacobien est donc une hypothèse raisonnable, mais si c'est le cas, elle est si subtile que je ne peux pas la voir. EDIT: J'ai modifié jacobian () pour utiliser les différences avant au lieu des différences centrales, et maintenant j'observe le développement attendu de l'erreur. Cependant, newton () ne parvient toujours pas à converger. En observant dx dans l'itération de Newton, je vois qu'il ne fait que croître, il n'y a même pas de fluctuation: il double presque (facteur 1.9) à chaque pas, le facteur augmentant progressivement.
Ascher et Petzold mentionnent que les approximations de différence pour le jacobien ne fonctionnent pas toujours bien. Un jacobien approximatif avec des différences finies peut-il provoquer une instabilité dans la méthode de Newton? Ou la cause est-elle ailleurs? Sinon, comment pourrais-je aborder ce problème?
Réponses:
Plus un long commentaire qu'autre chose:
Regardez le code pour l'approximation du quotient de différence de SUNDIALS pour avoir une meilleure idée de ce que vous devez faire dans une implémentation. Ascher et Petzold est un bon livre pour commencer, mais SUNDIALS est en fait utilisé dans le travail de production et a donc été mieux testé. (De plus, SUNDIALS est lié à DASPK, sur lequel Petzold a travaillé.)
Empiriquement, les Jacobiens approximatifs peuvent provoquer des échecs de convergence dans la méthode de Newton. Je ne sais pas si je les qualifierais d '«instabilité»; dans certains cas, il n'est tout simplement pas possible d'atteindre les tolérances d'erreur souhaitées dans les critères de terminaison. Dans d'autres cas, cela pourrait se manifester par une instabilité. Je suis presque certain qu'il y a un résultat plus quantitatif sur ce phénomène dans le livre des méthodes numériques de Higham, ou dans la discussion de Hairer et Wanner sur les méthodes W.
Cela dépend où vous pensez que l'erreur pourrait être. Si vous êtes extrêmement confiant dans votre implémentation d'Euler en arrière, je ne commencerais pas par là. L'expérience m'a rendu paranoïaque dans mes implémentations de méthodes numériques, donc si c'était moi, je commencerais par coder quelques problèmes de test vraiment basiques (quelques problèmes linéaires non rigides et rigides, l'équation de la chaleur avec une approximation centrée par différence finie, des choses comme ça) et j'utiliserais la méthode des solutions fabriquées pour m'assurer que je sais ce que sera la solution et ce que je devrais comparer.
Cependant, vous en avez déjà fait une partie:
Ce serait la prochaine chose que je testerais: utiliser un jacobien analytique. Après cela, vous pourriez également regarder les valeurs propres extrêmes de votre Jacobian à différence finie au cas où vous vous trouveriez dans la région instable d'Euler arriéré. L'examen des valeurs propres extrêmes de votre jacobien analytique comme base de comparaison pourrait vous donner un aperçu. En supposant que tous les vérifient, le problème est probablement dans la solution de Newton.
la source