Numerisch stabile explizite Lösung eines kleinen linearen Systems

Ich habe ein inhomogenes lineares System

wobei $A$ eine reelle $n\times n$ Matrix mit $n\leq 4$ . $A$ wird garantiert, dass der Nullraum von A eine Dimension von Null hat, sodass die Gleichung eine eindeutige Umkehrung von $x=A^{-1} b$ . Da das Ergebnis auf der rechten Seite einer ODE eingeht, die ich mit einer adaptiven Methode lösen möchte, ist es wichtig, dass die Lösung in Bezug auf kleine Variationen der Elemente von $A$ und glatt ist $b$. Aufgrund dieser Anforderung und der geringen Dimensionalität dachte ich, explizite Formeln für A - 1 b zu implementieren$A^{-1} b$. Die Elemente können genau Null sein oder sehr unterschiedliche Werte annehmen. Meine Frage ist, ob dies für Sie sinnvoll ist und ob dafür stabile Ausdrücke bekannt sind. Ich codiere in C für x86-Systeme.

Bevor ich explizite Formeln implementiere, würde ich mir die Frage stellen: "Lohnt es sich?":

• Lohnt es sich, diese expliziten Formeln zu schreiben, zu debuggen und zu validieren, während Sie leicht auf BLAS + LAPACK verlinken können, das die klassische Gaußsche Eliminierung verwendet?
• Erwarten Sie Stabilität? Ich glaube nicht, dass das Programmieren expliziter Formeln (wie Cramers Regel) im Gegenteil zu einer besseren Stabilität führt.
• Erwarten Sie, Geschwindigkeit zu gewinnen? Haben Sie bereits Ihr gesamtes Programm profiliert? Wie viel Zeit wird für die Lösung dieser 4x4-Systeme aufgewendet?
• Wie hoch ist die Wahrscheinlichkeit, dass Sie in einem Jahr Ihr Modell verbessern und 5 statt 4 Gleichungen benötigen?

Mein Rat: Verwenden Sie zuerst die BLAS / LAPACK-Kombination, prüfen Sie, ob sie funktioniert, profilieren Sie das gesamte Programm, bitten Sie einen Schüler, explizite Formeln zu implementieren (sorry, hier sarkastisch zu sein) und vergleichen Sie Geschwindigkeit und Robustheit.

$\mathcal{O}(n^{3})$

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

Um sicher zu gehen, ist es wahrscheinlich am besten sicherzustellen, dass keinen numerischen Rangmangel aufweist (dh keine kleinen Singularwerte hat).$A$

Das Problem mit der Cramer-Regel ist, dass ihre Stabilitätseigenschaften unbekannt sind, mit Ausnahme von (das vorwärts stabil, aber nicht rückwärts stabil ist). (Siehe Genauigkeit und Stabilität numerischer Algorithmen , 2. Auflage, von N. Higham.) Es wird nicht als zuverlässiger Algorithmus angesehen. Die Gaußsche Eliminierung mit partiellem Schwenken (GEPP) wird bevorzugt.$n = 2$

Ich würde erwarten, dass das Problem bei der Verwendung von BLAS + LAPACK zur Ausführung von GEPP in einer ODE-Lösung eine endliche Differenzierung ist, die in einer impliziten ODE-Methode verwendet wird. Ich weiß, dass Menschen lineare Programme im Rahmen einer Auswertung auf der rechten Seite gelöst haben, und weil sie dies naiv getan haben (nur die lineare Programmlösung in die rechte Seite gesteckt und einen Simplex-Algorithmus aufgerufen), haben sie die Genauigkeit ihrer Programme stark reduziert berechnete Lösung und wesentlich länger die Zeit, die es dauerte, um das Problem zu lösen. Ein Labmate von mir fand heraus, wie man solche Probleme viel effizienter und genauer lösen kann. Ich muss nachsehen, ob seine Veröffentlichung bereits veröffentlicht wurde. Möglicherweise haben Sie ein ähnliches Problem, unabhängig davon, ob Sie sich für die Verwendung von GEPP oder Cramer's Rule entscheiden.

Wenn Sie auf irgendeine Weise eine analytische Jacobi-Matrix für Ihr Problem berechnen können, möchten Sie dies möglicherweise tun, um sich einige numerische Kopfschmerzen zu ersparen. Es ist billiger zu bewerten und wahrscheinlich genauer als eine Näherung mit endlichen Differenzen. Ausdrücke für die Ableitung der Matrixinverse finden Sie hier, wenn Sie sie benötigen. Die Auswertung der Ableitung der Matrixinverse sieht so aus, als würde es mindestens zwei oder drei lineare Systemlösungen erfordern, aber sie würden alle dieselbe Matrix und unterschiedliche rechte Seiten haben, so dass sie nicht wesentlich teurer wäre als ein einzelnes lineares System lösen.

Und wenn es eine Möglichkeit gibt, Ihre berechnete Lösung mit einer Lösung mit bekannten Parameterwerten zu vergleichen, würde ich das tun, damit Sie diagnostizieren können, ob Sie auf eine dieser numerischen Fallstricke gestoßen sind.

Ich bin mir nicht sicher, ob das helfen kann, aber ich denke nur, wenn Sie über eine stabile Lösung sprechen, sprechen Sie über Approximationsmethoden. Wenn Sie Dinge explizit berechnen, hat Stabilität keinen Sinn. Das heißt, Sie müssen eine Näherungslösung akzeptieren, wenn Sie an Stabilität gewinnen möchten.