Solution explicite numériquement stable d'un petit système linéaire

J'ai un système linéaire inhomogène

où $A$ est une matrice $n\times n$ réelle avec $n\leq 4$ . L'espace nul de $A$ est garanti d'être de dimension nulle, l'équation a donc un inverse unique $x=A^{-1} b$ . Puisque le résultat entre dans le côté droit d'une ODE, que j'ai l'intention de résoudre en utilisant une méthode adaptative, il est important que la solution soit lisse en ce qui concerne les petites variations des éléments de $A$ et $b$ . En raison de cette exigence et de la petite dimensionnalité, j'ai pensé à mettre en œuvre des formules explicites pour $A^{-1} b$. Les éléments peuvent être exactement nuls ou prendre des valeurs très différentes. Ma question est de savoir si cela a du sens pour vous et s'il existe des expressions stables connues pour cela. Je code en C pour les systèmes x86.

Avant d'implémenter des formules explicites, je me pose la question: "ça vaut le coup?":

• Vaut-il la peine de passer du temps à écrire, déboguer et valider ces formules explicites alors que vous pourriez facilement vous lier à BLAS + LAPACK qui utilise l'élimination gaussienne classique?
• Vous attendez-vous à gagner en stabilité? Je ne pense pas que la programmation de formules explicites (comme la règle de Cramer) vous donnera au contraire une meilleure stabilité.
• Vous attendez-vous à gagner en vitesse? Avez-vous déjà profilé l'ensemble de votre programme? Quelle fraction de temps est consacrée à la résolution de ces systèmes 4x4?
• Quelle est la probabilité qu'en un an, vous amélioriez votre modèle et que vous ayez besoin de 5 équations au lieu de 4?

Mon conseil: utilisez d'abord la combinaison BLAS / LAPACK, voyez si cela fonctionne, profilez l'ensemble du programme, demandez à un étudiant d'implémenter des formules explicites (désolé, étant sarcastique ici) et faites une comparaison sur la vitesse et la robustesse.

Le seul résultat inverse explicite que je connaisse est la règle de Cramer , qui s'est récemment révélée calculable en temps (comme l'élimination gaussienne; cependant, incertain de la constante devant le facteur principal).$\mathcal{O}(n^{3})$

$A$$A$$\det(A) \neq 0$$x$$b$$x$$A$

Pour être sûr, il est probablement préférable de s'assurer que n'est pas non plus numériquement déficient en rang (c'est-à-dire qu'il n'a pas de petites valeurs singulières).$A$

Le problème avec la règle de Cramer est que ses propriétés de stabilité sont inconnues sauf pour (qui est stable en avant, mais pas en arrière). (Voir Exactitude et stabilité des algorithmes numériques , 2e édition, par N. Higham.) Il n'est pas considéré comme un algorithme fiable; L'élimination gaussienne avec pivot partiel (GEPP) est favorisée.$n = 2$

Je m'attendrais à ce que le problème avec BLAS + LAPACK pour effectuer GEPP dans une résolution ODE soit une différence finie utilisée dans une méthode ODE implicite. Je sais que les gens ont résolu des programmes linéaires dans le cadre d'une évaluation à droite, et parce qu'ils l'ont fait naïvement (il suffit de brancher le programme linéaire à résoudre sur le côté droit, en appelant un algorithme simplex), ils ont considérablement réduit la précision de leur solution calculée et a considérablement augmenté le temps nécessaire pour résoudre le problème. Un de mes collègues a trouvé un moyen de résoudre ces problèmes de manière beaucoup plus efficace et précise; Je vais devoir regarder pour voir si sa publication a encore été publiée. Vous pouvez avoir un problème similaire, que vous choisissiez d'utiliser GEPP ou la règle de Cramer.

S'il existe un moyen de calculer une matrice jacobienne analytique pour votre problème, vous souhaiterez peut-être le faire pour vous éviter quelques maux de tête numériques. Il sera moins coûteux à évaluer et probablement plus précis qu'une approximation par différence finie. Les expressions pour la dérivée de l'inverse de la matrice peuvent être trouvées ici si vous en avez besoin. L'évaluation de la dérivée de l'inverse de la matrice semble nécessiter au moins deux ou trois résolutions de système linéaire, mais elles seraient toutes avec la même matrice et différents côtés à droite, donc ce ne serait pas beaucoup plus cher qu'un système linéaire unique résoudre.

Et s'il existe un moyen de comparer votre solution calculée à une solution avec des valeurs de paramètres connues, je le ferais, afin que vous puissiez diagnostiquer si vous avez rencontré l'un de ces pièges numériques.

Je ne suis pas sûr que cela puisse aider, mais je pense simplement que lorsque vous parlez de solution stable, vous parlez de méthodes d'approximation. Lorsque vous calculez des choses explicitement, la stabilité n'a pas de sens. Cela signifie que vous devez accepter une solution approximative si vous voulez gagner en stabilité.