Reading Time: 12 minutes

Lorsque vous discréditez une équation aux dérivées partielles non linéaires, vous obtenez un système algébrique non linéaire. Pour le résoudre, il faut choisir la bonne méthode de solveur. La méthode de Newton converge quadratiquement mais exige un plein Jacobien. La méthode de Broyden construit le Jacobien moins cher mais peut déstabiliser. L’approche de Newton-Krylov sans Jacobien (JFNK) évite entièrement le Jacobien, ne s’appuyant que sur des produits à vecteur jacobien avec un préconditionnement basé sur la physique pour garder la convergence tractable.

Le choix n’est pas théorique. Il détermine si votre simulation converge du tout, combien de mémoire elle consomme et si elle se termine en quelques minutes ou en jours.

Lien vers l’article existant : Méthodes d’intégration du temps pour les solveurs PDE Couvre le paysage plus large des schémas explicites par rapport aux schémas implicites et de la division des opérateurs. Lien vers l’article existant : méthodes implicites ou explicites Discute des régions de stabilité, des critères de précision et lorsque des schémas implicites deviennent nécessaires. Cet article se concentre spécifiquement sur les solveurs non linéaires qui se trouvent dans chaque schéma implicite.

Lorsque les méthodes implicites exigent un solveur non linéaire

L’intégration implicite du temps nécessite la résolution d’un système non linéaire à chaque pas de temps. Considérez la discrétisation d’Euler rétrograde d’un PDE générique :

F(uⁿ⁺¹) = 0

L’inconnu uⁿ⁺¹ apparaît à l’intérieur de l’opérateur F par des termes de diffusion, une cinétique de réaction ou des équations multi-physiques couplées. Après discrétisation spatiale (volume fini, élément fini ou différence finie), le système devient :

F(x) = 0, where x ∈ ℝⁿ

La dimension n est le nombre de degrés de liberté – généralement le produit des cellules de la grille, des champs physiques et des itérations de pas de temps. Dans une simulation en champ de phase 3D avec 500 ³ de cellules et deux champs couplés, n peut facilement dépasser 10 ⁸.

Le défi est que F est non linéaire. Vous ne pouvez pas écrire x = F⁻¹(0). Au lieu de cela, vous avez besoin d’un algorithme itératif qui construit une séquence convergeant vers la racine.

La méthode de Newton : la ligne de base

La méthode de Newton est la référence par rapport à laquelle tous les autres solveurs non linéaires sont mesurés. Pour le système F(x) = 0, il construit la matrice jacobienne :

J = ∂F/∂x

et itère :

xₖ₊₁ = xₖ − J⁻¹·F(xₖ)

La propriété clé est convergence quadratique. Si le Jacobien est continu de Lipschitz et que la supposition initiale est dans le bassin de convergence, l’erreur satisfait :

‖xₖ₊₁ − x*‖ ≤ C · ‖xₖ − x*‖²

Cela signifie que le nombre de chiffres corrects double à chaque itération. Pour les systèmes PDE bien élevés, la méthode de Newton converge souvent en 3 à 5 itérations.

Jacobien complet : différence analytique ou finie

La construction du Jacobien est la partie la plus chère de la méthode de Newton. Les options sont :

Jacobien analytique. Vous dérivez ∂F/∂x à la main ou par la différenciation automatique. Dans des bibliothèques comme Fenics, la formulation variationnelle génère automatiquement le Jacobien par la différenciation symbolique. Le résultat est exact (jusqu’à une précision en virgule flottante).

Différence finie Jacobien. Vous perturbez chaque colonne de J par ε et réévaluez F :

Jᵢ ≈ (F(x + ε·eᵢ) − F(x)) / ε

Cela nécessite n des évaluations de fonctions supplémentaires par itération newton. Pour les grands n, cela domine le coût.

Bloquer Jacobian. Dans les systèmes couplés avec k champs, le Jacobien est structuré en bloc :

J = [ ∂F₁/∂x₁  ∂F₁/∂x₂  ... ]
    [ ∂F₂/∂x₁  ∂F₂/∂x₂  ... ]
    [          ...       ... ]

Si vous n’avez besoin que de blocs diagonaux ∂Fᵢ/∂xᵢ, vous pouvez les évaluer indépendamment – une technique appelée gel diagonal ou solveurs entièrement couplés ou séparés.

Convergence mondiale : recherche de lignes et régions de confiance

La méthode de Pure Newton est une méthode locale. Si la supposition initiale est loin de la solution, les itérations peuvent diverger même si la convergence quadratique est proche de la racine.

Recherche de ligne modifie l’étape Newton en la mettant à l’échelle :

xₖ₊₁ = xₖ − α·J⁻¹·F(xₖ)

où α ∈ (0, 1] est choisi pour réduire une fonction de mérite. L’approche classique minimise ‖F(xₖ − α·J⁻¹·F(xₖ))‖.

Les régions de confiance définissent un rayon Δ dans lequel le modèle quadratique est approuvé :

min q(s) = F(xₖ) + J·s + ½·sᵀ·H·s
subject to ‖s‖ ≤ Δ

Lorsque la réduction de la fonction de mérite est insuffisante, le rayon se rétrécit. Quand c’est suffisant, le rayon grandit.

Le package SNES de PETSC (solveur d’équations non linéaires SNU) fournit à la fois la mondialisation de recherche en ligne et de région de confiance. Le choix affecte le nombre total d’itérations plus que le coût par itération.

Newton pratique : quand ça marche et quand ça ne marche pas

La méthode de Newton excelle lorsque :

  • Le Jacobien est disponible sous une forme exacte (symbolique, AD ou bloc-analytique)
  • Les non-linéarités sont légères (diffusion, termes de réaction légers)
  • La mémoire est suffisante pour le Jacobien complet (ou sa factorisation clairsemée)

La méthode de Newton est difficile lorsque :

  • n dépasse 10 ⁶ et l’assemblage complet de Jacobien est prohibitif
  • La fonction F est opaque (code de simulation en boîte noire)
  • Chaque évaluation de fonction domine déjà le temps d’exécution

Pour les PDE à grande échelle, l’approche jacobienne complète est souvent trop chère. C’est là que les méthodes quasi-newton entrent.

Les méthodes quasi-Newton de Broyden

La méthode de Broyden appartient à la classe plus large des méthodes quasi-Newton. Au lieu de calculer J = ∂F/∂x, il crée une approximation Bₖ ≈ J par les mises à jour de rang 1 de la condition sécante :

Bₖ₊₁ · (xₖ₊₁ − xₖ) = F(xₖ₊₁) − F(xₖ)

La condition sécante est une équation pour n inconnue dans Bₖ₊₁. La formulation de Broyden choisit la mise à jour qui minimise le changement de norme Frobenius :

min ‖Bₖ₊₁ − Bₖ‖_F
subject to Bₖ₊₁·(xₖ₊₁ − xₖ) = F(xₖ₊₁) − F(xₖ)

bon contre mauvais broyden

La méthode d’origine de Broyden (parfois appelée « bon Broyden ») est mise à jour Bₖ à l’aide de l’équation sécante. La variante inverse (« Bad Broyden ») est mise à jour directement Bₖ⁻¹ :

Bₖ₊₁⁻¹ = Bₖ⁻¹ + (δy − Bₖ⁻¹·δx)·δxᵀ / (δxᵀ·δx)

où δx = xₖ₊₁ − xₖ et δy = F(xₖ₊₁) − F(xₖ).

La mise à jour inverse est ce que Python utilise scipy.optimize.broyden1. Cela évite d’inverser l’approximation jacobienne à chaque étape – une économie importante pour les gros problèmes.

L’instabilité du nombre de conditions

La limitation pratique la plus importante de la méthode de Broyden est que son approximation jacobienne peut devenir mal conditionnée. La mise à jour de rang 1 ne modifie qu’une seule direction, tandis que le reste de la matrice dérive. Au cours de nombreuses itérations, Bₖ accumule des informations de courbure obsolètes.

Le nombre de conditions κ(Bₖ) = ‖Bₖ‖ · ‖Bₖ⁻¹‖ peut croître sans borne. Un nombre de conditions excédant 10⁸ provoque le blocage des solveurs de Krylov comme les GMRES, même si le problème sous-jacent est bien conditionné.

C’est la raison concrète de préférer Newton ou JFNK à Broyden dans les simulations de production. Pour les problèmes légèrement non linéaires avec peu de degrés de liberté par champ (par exemple, diffusion 1D, régime permanent 2D), la méthode de Broyden peut être compétitive. Pour les systèmes 3D avec 10⁶+ degrés de liberté, la dérive de condition la rend peu fiable.

Broyden à mémoire limitée (L-BFGS)

Pour une optimisation à grande échelle, L-BFGS stocke les dernières paires de pas m (δxᵢ, δyᵢ) et construit implicitement l’inverse de Jacobien à travers une récursivité. Cela évite de stocker la matrice complète n × n.

Cependant, L-BFGS est conçu pour l’optimisation (minimation d’un objectif scalaire), et non la recherche de racine (résolution F(x) = 0). Pour la recherche de racines PDE, la mise à jour de Broyden à mémoire limitée utilisée dans les solveurs comme le solveur BFGS de PETSC (et non L-BFGS) est plus courante.

Newton-Krylov sans Jacobian : la norme à grande échelle moderne

La méthode Newton-Krylov (JFNK) sans Jacobien aborde le goulot d’étranglement de la mémoire de Jacobien Newton en assemblant jamais le Jacobien. Au lieu de cela, il utilise une méthode de sous-espace Krylov (généralement GMRES) pour résoudre le système Newton :

J · Δx = −F(x)

GMRES ne nécessite que les produits Matrix-Vector J·v. Ceux-ci sont calculés via le produit Jacobian-Vector :

J · v ≈ (F(x + ε·v) − F(x)) / ε

Aucun Jacobien explicite n’est assemblé. Seuls F et le vecteur de perturbation v sont nécessaires.

Pourquoi le préconditionnement est obligatoire

JFNK avec un produit à vecteur jacobien nu est rarement pratique. Le système de Newton peut avoir un grand nombre de conditions, en particulier lorsque le Jacobien a des valeurs propres très variables à partir de différents processus physiques.

Préconditionnement basé sur la physique (PBP) aborde ce problème en tenant compte des opérateurs connus. Envisagez un système de convection-diffusion-réaction :

J = J_convection + J_diffusion + J_reaction

PBP remplace le préconditionneur par :

M = J_diffusion

L’opérateur de diffusion domine généralement le spectre et est souvent tridiagonal ou tridiagonal par champ. Son inverse peut être calculé efficacement par la factorisation directe clairsemée ou l’ILU, tandis que les termes de convection et de réaction sont traités explicitement.

Le système préconditionné devient :

M⁻¹ · J · Δx = −M⁻¹ · F(x)

Le nombre de conditions de M⁻¹·J est considérablement réduit car l’opérateur de diffusion – la source du pire conditionnement – a été pris en compte dans M.

Sans préconditionnement, JFNK peut nécessiter des centaines ou des milliers d’itérations GMRES. Avec le PBP, il converge souvent en 10 à 50 itérations, correspondant ou dépassant le Newton Jacobien.

Le modèle Newton-Krylov en pratique

  1. Construire le résidu newton r = F(xₖ)
  2. Résoudre J · Δx = −r Utiliser GMRES
  3. Les produits Jacobian-Vector sont calculés par des perturbations de différences finies de F
  4. Preconditionneur M est appliqué à chaque itération GMRES
  5. Mise à jour xₖ₊₁ = xₖ + Δx
  6. Vérifier la convergence sur ‖r‖

Ce modèle apparaît dans tous les principaux cadres PDE : FENICS/Dolfin, Petsc/SNES, FIPY et NonLinearSolve.jl. La différence réside dans la façon dont le préconditionneur est construit et la façon dont le produit de vecteur jacobien est calculé.

Le modèle de polyalgorithme : comment les solveurs modernes se déploient réellement

Les solveurs de production s’appuient rarement sur une seule méthode. Au lieu de cela, ils utilisent un polyalgorithm qui s’adapte à la difficulté du problème :

1. Start with fast Broyden (cheap per iteration, no Jacobian)
2. If convergence stalls, fall back to full Newton (robust quadratic convergence)
3. If Newton also fails, fall back to TrustRegion (globally convergent)

Il s’agit du modèle utilisé par PETSC SNES, DifferentialEquations.jl de Sciml et le solveur non linéaire de Moose.

Eisenstat-Walker forçant

La stratégie Eisenstat-Walker (parfois appelée la « stratégie de force ») adapte la tolérance intérieure du solveur de Krylov par rapport au résidu de Newton :

‖J · Δx + r‖ ≤ ηₖ · ‖r‖

où ηₖ se détend comme r rétrécit. Les premières itérations de Newton utilisent une tolérance intérieure lâche (moins de pas de Krylov). Au fur et à mesure que le résidu de Newton diminue, la tolérance intérieure se resserre (plus de précision de Krylov).

Cela empêche les travaux perdus de résoudre le système linéaire avec une grande précision lorsque le Newton Iterate est encore loin de la racine.

Cadre de décision : quand choisir la méthode

Le tableau suivant résume les compromis pratiques :

Comparaison des méthodes

Critère Newton-Raphson Broyden (quasi-Newton) JFNK (Newton-Krylov)
taux de convergence Quadratique (racine proche) Superlinéaire (par itération), mais peu fiable pour les grands n Dépend du préconditionneur ; Super linéaire avec un bon PBP
Empreinte mémoire O(n²) pour Jacobien complet (mais le peu moins le réduit) O(n²) pour l’approximation inverse de Broyden ; O(m·n) pour les variantes à mémoire limitée O(n) — Uniquement les vecteurs de fonction et le stockage de Krylov
Coût par itération Élevé (assemblage jacobien + factorisation) Faible (mise à jour de rang 1 + évaluation de la fonction) Modéré (JVP + itération Krylov)
Exigence Jacobienne Jacobien complet (analytique, AD ou FD) Pas de Jacobien (approprié à partir des étapes) Pas de Jacobien explicite (JVP uniquement)
préconditionnement Facultatif (améliore la résolution linéaire à Newton) Non applicable (sans jacobie) Requis pour la convergence pratique (PBP fortement recommandé)
Meilleur cas d’utilisation Problèmes de petite à moyenne avec Jacobien exact disponible Problèmes 1D-2D légèrement non linéaires avec un budget limité Des problèmes 3D à grande échelle avec des millions de degrés de liberté

Choisir par taille de problème

Petit problèmes (n < ; 10 ⁴) : Newton complet avec jacobien analytique ou symbolique. Le coût est dominé par l’assemblage matriciel, de sorte que la convergence quadratique l’emporte sur le coût de configuration.

Problèmes moyens (10⁴ < n < 10⁶): Broyden peut être compétitif si la non-linéarité est légère. Newton avec une factorisation directe clairsemée est plus robuste mais coûte plus cher par itération.

De gros problèmes (n > 10⁶): JFNK avec un préconditionnement basé sur la physique est le choix standard. L’assemblage jacobien complet est prohibitif et la dérive de l’état de Broyden n’est pas fiable.

Choisir par non-linéarité

Non linéarité légère (termes de réaction faiblement dominés par la diffusion) : Newton ou Broyden fonctionnent bien. La convergence quadratique de Newton le rend attrayant pour les petits systèmes.

La non-linéarité forte (champ de phase couplé, réaction-diffusion avec une cinétique rigide) : Newton avec la mondialisation (recherche de ligne ou région de confiance) est essentielle. L’approximation de Broyden peut induire la direction de recherche en erreur.

Non-linéarité opaque (simulation en boîte noire, inconnue F): JFNK ou Broyden, car ils ne nécessitent pas de construction jacobienne explicite.

Choisir par des contraintes de mémoire

Mémoire limitée (single-nœud, < 64 Go) : JFNK ne nécessite que le stockage O(n) pour les vecteurs. Newton nécessite O(nz(j)) pour le stockage jacobien clairsemé plus O(nNZ(L)) pour la factorisation.

Mémoire complète (noeud de cluster, > 256 Go) : Newton complet devient réalisable pour des problèmes avec N jusqu’à 10 ⁷, à condition que le Jacobien soit clairsemé.

Exemples de code : à quoi ressemblent ces méthodes en pratique

Scipy : la méthode de Broyden

Scipy fournit broyden1 pour la recherche de racine avec la mise à jour quasi-Newton de Broyden :

from scipy.optimize import broyden1, root

def f(x):
    """System F(x) = 0"""
    return [x[0]**2 + x[1] - 1,
            x[0] + x[1]**2 - 2]

result = broyden1(f, [0.5, 0.5])
print(f"Solution: {result.x}")

La fonction broyden1 construit l’approximation jacobienne inverse grâce aux mises à jour de rang 1. Pour les petits systèmes, il évite l’assemblage complet de Jacobien. Pour les grands systèmes, la dérive de l’état peut entraîner une défaillance de la convergence.

L’interface root prend également en charge 'hybr' (Levenberg-Trust-région de Minpack) et 'lm' (Levenberg-Marquardt), qui peuvent tous deux mieux gérer les problèmes structurés des moindres carrés que Broyden.

Fenics : Newton avec Jacobien symbolique

Fenics génère symboliquement le Jacobien à partir de la forme variationnelle :

from fenics import *

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "P", 1)

u = TrialFunction(V)
v = TestFunction(V)

f = Constant(1.0)
a = dot(grad(u), grad(v)) * dx
L = f * v * dx

u_solution = Function(V)
solve(a == L, u_solution)

# Nonlinear: Newton iterate with automatic Jacobian
U = Function(V)  # current solution
u = TrialFunction(V)

F = dot(grad(U), grad(v)) * dx - f * v * dx  # residual
A = derivative(F, U, u)  # Jacobian (symbolic)

problem = NonlinearProblem(F, U, bcs)
solver = NewtonSolver(MPI.comm.world)
solver.parameters["linear_solver"] = "petsc"
solver.solve(problem)

L’appel derivative(F, U, u) génère le Jacobien exact grâce à la différenciation automatique. C’est la caractéristique la plus puissante de Fenics pour les problèmes non linéaires.

PETSC/SNES : Newton-Krylov avec préconditionnement

Le package SNES de PETSC met en œuvre le modèle complet de Newton-Krylov :

from pysns import SNES

# Define residual function
def residual(x):
    # Return F(x) — the nonlinear residual
    return F_of_x(x)

# Create SNES solver
snes = SNESCreate()
snesSetFunction(snes, x, residual)

# Set Newton-Krylov with preconditioning
SNESSetType(snes, SNESNEWTONKRYLOV)
SNESSetKrylovDimension(snes, 50)  # GMRES max iterations

# Configure preconditioner (physics-based)
ksp = SNESGetKSP(snes)
KSPSetType(ksp, KSPPRECONDEL)  # ILU preconditioner
KSPSetPreconditioner(ksp, PCILU)

# Solve
SNESSetUp(snes)
SNESSolve(snes)

Le type Newton-Krylov bascule automatiquement entre les variantes complètes de Newton (assemblée Jacobienne), de Newton sans Jacobian et de Broyden en fonction de la configuration SNESSetUser Jacobian. La configuration du préconditionneur contrôle le taux de convergence de Krylov.

NonLinearSolve.jl : le motif polyalgorithme

L’écosystème des différenciations de Julia utilise un polyalgorithme qui enchaîne les méthodes :

using DifferentialEquations

# Define the nonlinear system
f! = (residual, x) -> begin
    residual[1] = x[1]^2 + x[2] - 1
    residual[2] = x[1] + x[2]^2 - 2
end

# Newton-Krylov with Eisenstein-Walker forcing
prob = NonlinearProblem(f!, [0.5, 0.5])
sol = solve(prob, NewtonRaphson(); abstol=1e-8, reltol=1e-6)

# For large-scale systems: Jacobian-free Newton-Krylov
sol = solve(prob, JFNK(); abstol=1e-8)

L’algorithme JFNK calcule les produits à vecteur jacobien par des différences finies et délègue à un solveur de sous-espace de Krylov. Pour les systèmes rigides, le polyalgorithme revient à TrustRegion si JFNK décroche.

Performances empiriques : JFNK vs Opérateur-Split

Une référence concrète de la littérature démontre l’avantage pratique de JFNK. Dans les calculs de transfert radiatif non-LTE (NLTE), une étude de 2024 en astronomie & L’astrophysique (A&a) a comparé Newton-Krylov sans Jacobiens à des approches de fractionnement d’opérateur. JFNK a convergé approximativement 2× que le fractionnement de l’opérateur pour la même cible de précision, tout en utilisant une mémoire comparable. Le résultat met en évidence que les solveurs non linéaires monolithiques peuvent surpasser les approches partitionnées lorsque la non-linéarité est forte et que le préconditionneur est bien adapté à la physique.

Résumé : Quelle méthode recommandons-nous ?

La recommandation dépend de l’échelle et de la structure de votre problème :

  • Petits problèmes avec Jacobien exact disponible : Newton-Raphson. La convergence quadratique et le bassin de convergence robuste l’emportent sur le coût jacobien.
  • Systèmes 1D-2D non linéaires légèrement non linéaires : La méthode de Broyden est rentable. Acceptez le risque de dérive des conditions pour les petits et moyens problèmes où chaque itération est bon marché.
  • Systèmes 3D à grande échelle avec des millions de degrés de liberté : JFNK avec un préconditionnement basé sur la physique. C’est la seule option viable lorsque l’assemblage complet de Jacobien est prohibitif et que la dérive de l’état de Broyden est inacceptable.
  • Codes de simulation en boîte noire : JFNK ou Broyden. Les deux évitent la construction jacobienne explicite.
  • Exigences de simulation critiques : Newton avec la mondialisation de recherche en ligne. La convergence quadratique et la stabilité mondiale des méthodes de la région de confiance valent le coût supplémentaire.

Pour la plupart des workflows PDE de production, le modèle de polyalgorithme – Fast Broyden repli vers Newton Fallback to TrustRegion – offre le meilleur équilibre entre vitesse, robustesse et fiabilité.

Guides connexes

Références

Vous voulez de l’aide pour mettre en place un solveur non linéaire pour votre recherche ?

Choisir la bonne méthode de solveur non linéaire peut faire la différence entre une simulation qui se termine en heures et une autre qui ne converge jamais. Si vous êtes aux prises avec la divergence de Newton Iteration, Broyden Condition Drift ou JFNK Preconditioner Tuning, notre équipe peut vous aider.

Nous nous spécialisons dans la création de cadres de vérification pour les codes Python et Julia scientifiques, y compris les workflows de solveurs non linéaires pour les simulations PDE. Contactez-vous via notre système de suivi des problèmes pour discuter des besoins de votre projet.