Löse eine Matrixgleichung nach Jacobis Methode (überarbeitet)

11

Mathematischer Hintergrund

Sei A eine N mal N Matrix von reellen Zahlen, ein ba Vektor von N reellen Zahlen und xa Vektor N unbekannte reelle Zahlen. Eine Matrixgleichung lautet Ax = b.

Jacobis Methode ist wie folgt: Zerlegen Sie A = D + R, wobei D die Matrix der Diagonalen ist und R die verbleibenden Einträge sind.

Wenn Sie eine anfängliche Schätzlösung x0 erstellen, ist eine verbesserte Lösung x1 = invers (D) * (b - Rx), wobei alle Multiplikationen Matrix-Vektor-Multiplikation sind und invers (D) die inverse Matrix ist.


Problemspezifikation

  • Eingabe : Ihr gesamtes Programm sollte die folgenden Daten als Eingabe akzeptieren: die Matrix A, den Vektor b, eine anfängliche Schätzung x0 und eine 'Fehlernummer' e.
  • Ausgabe : Das Programm muss die Mindestanzahl von Iterationen so ausgeben, dass sich die neueste Lösung von der tatsächlichen Lösung unterscheidet, höchstens von e. Dies bedeutet, dass sich jede Komponente der Vektoren in absoluter Größe um höchstens e unterscheidet. Sie müssen Jacobis Methode für Iterationen verwenden.

Wie die Daten eingegeben werden, liegt bei Ihnen . Es könnte Ihre eigene Syntax in einer Befehlszeile sein, Sie könnten Eingaben aus einer Datei übernehmen, was auch immer Sie wählen.

Wie die Daten ausgegeben werden, liegt bei Ihnen . Es kann in eine Datei geschrieben werden, die in der Befehlszeile angezeigt wird und als ASCII-Grafik geschrieben ist, solange es von einem Menschen gelesen werden kann.

Weitere Details

Sie erhalten nicht die wahre Lösung: Wie Sie die wahre Lösung berechnen, liegt ganz bei Ihnen. Sie können es beispielsweise nach der Cramer-Regel lösen oder eine Inverse direkt berechnen. Was zählt, ist, dass Sie eine echte Lösung haben, um mit Iterationen vergleichen zu können.

Präzision ist ein Problem; Die "genauen Lösungen" einiger Leute zum Vergleich können unterschiedlich sein. Für die Zwecke dieses Code-Golfs muss die genaue Lösung 10 Dezimalstellen entsprechen.

Um ganz klar zu sein: Wenn sogar eine Komponente Ihrer aktuellen Iterationslösung die entsprechende Komponente in der wahren Lösung um e überschreitet, müssen Sie die Iteration fortsetzen.

Die Obergrenze für N hängt davon ab, welche Hardware Sie verwenden und wie viel Zeit Sie für die Ausführung des Programms aufwenden möchten. Für die Zwecke dieses Code-Golfs wird maximal N = 50 angenommen.

Voraussetzungen

Wenn Ihr Programm aufgerufen wird, können Sie jederzeit davon ausgehen, dass Folgendes gilt:

  • N> 1 und N <51, dh Sie erhalten niemals eine Skalargleichung, immer eine Matrixgleichung.
  • Alle Eingaben liegen über dem Feld der reellen Zahlen und werden niemals komplex sein.
  • Die Matrix A ist immer so, dass die Methode zur wahren Lösung konvergiert, so dass Sie immer eine Anzahl von Iterationen finden können, um den Fehler (wie oben definiert) unten oder gleich e zu minimieren.
  • A ist niemals die Nullmatrix oder die Identitätsmatrix, dh es gibt eine Lösung.

Testfälle

A = ((9, -2), (1, 3)), b = (3,4), x0 = (1,1), e = 0.04

Die wahre Lösung ist (0,586, 1,138). Die erste Iteration ergibt x1 = (5/9, 1), was sich um mehr als 0,04 von der tatsächlichen Lösung um mindestens eine Komponente unterscheidet. Bei einer weiteren Iteration ergibt sich x2 = (0,555, 1,148), was sich um weniger als 0,04 von (0,586, 1,138) unterscheidet. Somit ist die Ausgabe

2

A = ((2, 3), (1, 4)), b = (2, -1), x0 = (2.7, -0.7), e = 1.0

In diesem Fall ist die wahre Lösung (2.2, -0.8) und die anfängliche Schätzung x0 hat bereits einen Fehler von weniger als e = 1.0, daher geben wir 0 aus. Das heißt, wenn Sie keine Iteration durchführen müssen, geben Sie einfach aus

0

Einreichungsbewertung

Dies ist Code Golf, wobei alle Standardlücken hiermit nicht zugelassen werden. Das kürzeste korrekte vollständige Programm (oder die Funktion), dh die niedrigste Anzahl von Bytes, gewinnt. Es wird davon abgeraten , Dinge wie Mathematica zu verwenden, die viele der erforderlichen Schritte in einer Funktion zusammenfassen, aber jede gewünschte Sprache verwenden.

user1997744
quelle
2
Sie sollten wirklich warten, um mehr Feedback zu erhalten, insbesondere angesichts des kürzlich geschlossenen Beitrags. PPCG-Herausforderungen haben normalerweise eine gemeinsame Struktur in den Spezifikationen, was normalerweise dazu beiträgt, dass sie leicht zu verstehen sind und nicht lästig und mehrdeutig. Versuchen Sie, einige der einigermaßen positiv bewerteten Herausforderungen zu betrachten und das Format nachzuahmen.
Uriel
@Uriel Ich erkenne dies, aber ich habe das Gefühl, dass meine Spezifikation vollständig ist, und das Format passt zwar nicht genau zu den meisten Fragen, kann jedoch linear gelesen werden und führt den Leser entsprechend. Das Format sollte auch den Inhalt des Problems selbst berücksichtigen.
user1997744
3
"Das kürzeste richtige vollständige Programm " klingt so, als würden Sie nur Programme und keine Funktionen zulassen. Ich würde "/ function" hinzufügen.
Adám
2
+1 Formatierung macht oder bricht die Fähigkeit meines Gehirns, sich auf eine Frage zu konzentrieren
Stephen
1
@ user1997744 Ja, macht Sinn. Ich glaube, die Standardeinstellung ist, dass jeder andere Code, wie andere Funktionen oder Python-Importe, zulässig ist, aber auch in der Anzahl der Bytes enthalten ist.
Stephen

Antworten:

4

APL (Dyalog) , 78 68 65 49 Bytes

Genau die Art des Problems, für das APL erstellt wurde.

-3 danke an Erik den Outgolfer . -11 danke an ngn .

Anonyme Infix-Funktion. Nimmt A als linkes Argument und x als rechtes Argument. Druckt das Ergebnis in STDOUT als vertikal unär mit 1Tally-Markierungen, gefolgt von 0Interpunktion. Dies bedeutet, dass sogar ein 0-Ergebnis sichtbar ist, das kein 1s vor dem ist 0.

{⎕←∨/e<|⍵-b⌹⊃A b e←⍺:⍺∇D+.×b-⍵+.×⍨A-⌹D←⌹A×=/¨⍳⍴A}

Probieren Sie es online aus!

Erklärung in Lesereihenfolge

Beachten Sie, wie der Code der Problemspezifikation sehr ähnlich lautet:

{…  Drucken Sie } auf dem gegebenen A, b und e und dem gegebenen x,  ob die Aussage, dass  e kleiner als  der Absolutwert von x minus  b ist, geteilt durch  den ersten von A, b und e (dh A),  die das linke Argument sind,  und wenn ja,  rekursieren Sie auf  D Matrixzeiten  b minus  x, Matrix multipliziert mit  A minus  der Umkehrung von D (die verbleibenden Einträge),  wobei D  A ist, wo  es gleiche  Koordinaten für  die Form gibt von A (dh die Diagonale)
⎕←
∨/
e<
|⍵-
b⌹
⊃A b e
←⍺
:
  ⍺∇
  D+.×
  b-
  ⍵+.×⍨
  A-
  ⌹D
  
  
  =/¨
  
  ⍴A

Schritt für Schritt Erklärung

Die tatsächliche Ausführungsreihenfolge von rechts nach links:

{} Anonyme Funktion, bei der A be und ⍵ x ist:
A b c←⍺ linkes Argument in A, b und e aufteilen Wählen Sie
 die erste (A)
b⌹ Matrixteilung mit b (gibt den wahren Wert von x an)
⍵- Differenzen zwischen aktuellen Werten von x und diesen  akzeptablen
| absoluten Werten
e<Fehler weniger als diese?
∨/ wahr für irgendeinen? (Lit. ODER Verkleinerung)
⎕← Drucken Sie diesen Booleschen Wert auf STDOUT
: und wenn ja:
  ⍴A Form einer
   Matrix dieser Form, bei der jede Zelle ihre eigenen Koordinaten
  =/¨ für jede Zelle hat. Sind die vertikalen und horizontalen Koordinaten gleich? (diagonal)
   multiplizieren Sie die Zellen von A mit dem
   inversen (diagonalen) Matrix- Inversspeicher
  D← in D (für D iagonal)
   invers (zurück zum Normalen)
  A- subtrahieren von A-
  ⍵+.×⍨ Matrix multiplizieren (dasselbe wie Punktprodukt, daher das .), dass mit x
  b- das von b-
  D+.× Matrixprodukt von D subtrahiert wird und dass
  ⍺∇ diese Funktion mit gegebenem A be angewendet wird und dass als neuer Wert von x

Adam
quelle
Die Ausgabe sollte die Anzahl der Iterationen sein, die für eine Genauigkeit von erforderlich sind e.
Zgarb
-1: Dies ist keine Lösung. Sie benötigen x0, da der springende Punkt darin besteht, zu wissen, wie viele Schritte erforderlich sind, um von einem bestimmten Startpunkt aus zu einer gewünschten Genauigkeit zu gelangen.
user1997744
@ user1997744 Ah, ich habe das Problem falsch verstanden. Es tut uns leid.
Adám
@ user1997744 Besser?
Adám
1
@ user1997744 Nicht arithmetische Operation, nur die Fähigkeit zu lesen einstellig , wo in der Tat 0 nichts .
Adám
1

Python 3 , 132 Bytes

f=lambda A,b,x,e:e<l.norm(x-dot(l.inv(A),b))and 1+f(A,b,dot(l.inv(d(d(A))),b-dot(A-d(d(A)),x)),e)
from numpy import*
l=linalg
d=diag

Probieren Sie es online aus!

Verwendet eine rekursive Lösung.

notjagan
quelle
@ Adám Ich bin mir nicht sicher, ob ich das wirklich verstehe. Ich habe es so interpretiert, fdass der Name nicht im Codeblock enthalten ist, den ich jetzt behoben habe. Wenn es sich jedoch um ein völlig anderes Problem handelt, kann es dennoch ein Problem sein.
Notjagan
@ Adám Diese Antwort scheint zu bestätigen, was ich derzeit habe; Es ist eine Funktion mit Hilfscode, die nach ihrer Definition als Einheit arbeiten kann.
Notjagan
Ah, ok. Egal Dann. Ich kenne Python nicht, also war ich nur neugierig. Gut gemacht!
Adám
Ist das Stoppkriterium nicht "Dies bedeutet, dass sich jede Komponente der Vektoren in absoluter Größe um höchstens e unterscheidet"? Grundsätzlich die Max-Norm, nicht die L2-Norm.
NikoNyrh
@NikoNyrh behoben.
Notjagan
1

R 138 Bytes

function(A,x,b,e){R=A-(D=diag(diag(A)))
g=solve(A,b)
if(norm(t(x-g),"M")<e)T=0
while(norm((y=solve(D)%*%(b-R%*%x))-g,"M")>e){T=T+1;x=y}
T}

Probieren Sie es online aus!

Vielen Dank an NikoNyrh für die Behebung eines Fehlers

Es ist auch erwähnenswert, dass es ein R-Paket gibt Rlinsolve, das eine lsolve.jacobiFunktion enthält , die eine Liste mit x(der Lösung) und iter(der Anzahl der erforderlichen Iterationen) zurückgibt, aber ich bin nicht sicher, ob es die richtigen Berechnungen durchführt.

Giuseppe
quelle
Ist das Stoppkriterium nicht "Dies bedeutet, dass sich jede Komponente der Vektoren in absoluter Größe um höchstens e unterscheidet"? Grundsätzlich die Max-Norm, nicht die L2-Norm.
NikoNyrh
@NikoNyrh du bist richtig! normZum Glück bietet die Funktion das auch für mich ohne zusätzliche Bytes.
Giuseppe
1

Clojure, 212 198 196 Bytes

#(let[E(range)I(iterate(fn[X](map(fn[r b d](/(- b(apply +(map * r X)))d))(map assoc % E(repeat 0))%2(map nth % E)))%3)](count(for[x I :while(not-every?(fn[e](<(- %4)e %4))(map -(nth I 1e9)x))]x)))

Ohne Matrixbibliothek implementiert, wiederholt es den Prozess 1e9 Mal, um die richtige Antwort zu erhalten. Dies würde bei zu schlecht konditionierten Eingängen nicht funktionieren, sollte aber in der Praxis gut funktionieren.

Weniger Golf, ich war ziemlich zufrieden mit den Ausdrücken von Rund D:) Die erste Eingabe %(A) muss ein Vektor sein, keine Liste, damit sie assocverwendet werden kann.

(def f #(let[R(map assoc %(range)(repeat 0))
             D(map nth %(range))
             I(iterate(fn[X](map(fn[r x b d](/(- b(apply +(map * r x)))d))R(repeat X)%2 D))%3)]
          (->> I
               (take-while (fn[x](not-every?(fn[e](<(- %4)e %4))(map -(nth I 1e9)x))))
               count)))
NikoNyrh
quelle