Résoudre une équation matricielle par la méthode de Jacobi (révisée)

11

Contexte mathématique

Soit A une matrice N par N de nombres réels, un vecteur ba de N nombres réels et xa vecteur N nombres réels inconnus. Une équation matricielle est Ax = b.

La méthode de Jacobi est la suivante: décomposer A = D + R, où D est la matrice des diagonales et R les entrées restantes.

si vous faites une solution initiale x0, une solution améliorée est x1 = inverse (D) * (b - Rx) où toutes les multiplications sont une multiplication matrice-vecteur et l'inverse (D) est l'inverse de la matrice.


Spécification du problème

  • Entrée : Votre programme complet doit accepter comme entrée les données suivantes: la matrice A, le vecteur b, une estimation initiale x0 et un numéro «d'erreur» e.
  • Sortie : Le programme doit sortir le nombre minimum d'itérations de telle sorte que la dernière solution diffère par la vraie solution, par au plus e. Cela signifie que chaque composante des vecteurs en amplitude absolue diffère d'au plus e. Vous devez utiliser la méthode de Jacobi pour les itérations.

La façon dont les données sont entrées est votre choix ; cela peut être votre propre syntaxe sur une ligne de commande, vous pouvez prendre des entrées à partir d'un fichier, quoi que vous choisissiez.

Le mode de sortie des données est votre choix ; il peut être écrit dans un fichier, affiché dans la ligne de commande, écrit en tant qu'art ASCII, n'importe quoi, tant qu'il est lisible par un humain.

Plus de détails

On ne vous donne pas la vraie solution: la façon dont vous calculez la vraie solution dépend entièrement de vous. Vous pouvez le résoudre par la règle de Cramer par exemple, ou en calculant directement un inverse. Ce qui compte, c'est que vous ayez une vraie solution pour pouvoir comparer aux itérations.

La précision est un problème; Les «solutions exactes» de comparaison de certaines personnes peuvent différer. Aux fins de ce code golf, la solution exacte doit être fidèle à 10 décimales.

Pour être absolument clair, si même un composant de votre solution d'itération actuelle dépasse son composant correspondant dans la vraie solution par e, vous devez continuer à itérer.

La limite supérieure de N varie en fonction du matériel que vous utilisez et du temps que vous êtes prêt à passer à exécuter le programme. Aux fins de ce code golf, supposons un maximum N = 50.

Conditions préalables

Lorsque votre programme est appelé, vous êtes libre de supposer que les éléments suivants sont toujours valables:

  • N> 1 et N <51, c'est-à-dire qu'on ne vous donnera jamais une équation scalaire, toujours une équation matricielle.
  • Toutes les entrées concernent le domaine des nombres réels et ne seront jamais complexes.
  • La matrice A est toujours telle que la méthode converge vers la vraie solution, de sorte que vous pouvez toujours trouver un certain nombre d'itérations pour minimiser l'erreur (telle que définie ci-dessus) ci-dessous ou égale à e.
  • A n'est jamais la matrice zéro ou la matrice d'identité, c'est-à-dire qu'il y a une solution.

Cas de test

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

La vraie solution est (0,586, 1,138). La première itération donne x1 = (5/9, 1), différant de plus de 0,04 de la vraie solution, d'au moins un composant. En prenant une autre itération, nous trouvons, x2 = (0,555, 1,148) qui diffère de moins de 0,04 de (0,586, 1,138). Ainsi, la sortie est

2

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

Dans ce cas, la vraie solution est (2.2, -0.8) et la supposition initiale x0 a déjà une erreur inférieure à e = 1.0, donc nous émettons 0. Autrement dit, chaque fois que vous n'avez pas besoin de faire une itération, vous émettez simplement

0

Évaluation de la soumission

Il s'agit du code golf, avec toutes les lacunes standard par la présente interdites. Le programme complet (ou la fonction) correct le plus court , c'est-à-dire le nombre d'octets le plus faible, gagne. Il est déconseillé d'utiliser des choses comme Mathematica qui regroupent un grand nombre des étapes nécessaires dans une seule fonction, mais utilisez la langue de votre choix.

user1997744
la source
2
Vous devriez vraiment attendre d'obtenir plus de commentaires à ce sujet, surtout compte tenu du récent message fermé. Les défis PPCG partagent généralement une structure commune dans les spécifications, ce qui contribue généralement à ce qu'ils soient faciles à comprendre, plutôt que fastidieux et ambigus. Essayez d'examiner certains des défis raisonnablement positifs et imitez le format.
Uriel
@Uriel Je m'en rends compte, mais je pense que j'ai été exhaustif dans mes spécifications, et le format, sans correspondre exactement à la majorité des questions, peut être lu de manière linéaire et guider le lecteur en conséquence. Le format doit également garder à l'esprit le contenu du problème lui-même.
user1997744
3
"Le programme complet correct le plus court " sonne comme si vous n'autorisez que les programmes et non les fonctions. J'ajouterais "/ fonction".
2017
2
La mise en forme +1 rend ou brise la capacité de mon cerveau à se concentrer sur une question
Stephen
1
@ user1997744 Oui, c'est logique. Je crois que la valeur par défaut est que tout autre code, comme les autres fonctions ou les importations python, est autorisé, mais également inclus dans le bytecount.
Stephen

Réponses:

4

APL (Dyalog) , 78 68 65 49 octets

Exactement le type de problème APL a été créé pour.

-3 merci à Erik l'Outgolfer . -11 grâce à ngn .

Fonction d'infixation anonyme. Prend A comme argument de gauche et x comme argument de droite. Les impressions résultent dans STDOUT sous forme d'unaire vertical en utilisant des 1marques de pointage, suivies d'une 0ponctuation. Cela signifie que même un résultat 0 peut être vu, étant aucun 1s avant le 0.

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

Essayez-le en ligne!

Explication dans l'ordre de lecture

Remarquez comment le code se lit de manière très similaire à la spécification du problème:

{} Sur les données A, b et e données et et sur le x donné,
⎕← indiquez
∨/ s'il existe une vérité dans l'énoncé selon lequel
e< e est inférieur à
|⍵- la valeur absolue de x moins
b⌹ b matrice divisée par
⊃A b e le premier de A, b et e (c'est-à-dire A)
←⍺ qui sont l'argument de gauche
: et si c'est le cas,
  ⍺∇ récapitulons sur
  D+.× D matrice fois
  b- b moins
  ⍵+.×⍨ x, matrice multipliée par
  A- A moins
  ⌹D l'inverse de D (les entrées restantes)
   où D est
   A où
  =/¨ il y a des
   coordonnées égales pour
  ⍴A la forme de A (c'est-à-dire la diagonale)

Explication étape par étape

L'ordre d'exécution réel de droite à gauche:

{} Fonction anonyme où est A be et ⍵ est x:
A b c←⍺ divise l'argument gauche en A, b et e
 sélectionne la première
b⌹ division de matrice (A) avec b (donne la vraie valeur de x)
⍵- différences entre les valeurs actuelles de x et les
| valeurs absolues
e< acceptables erreur moins que celles-ci?
∨/ vrai pour tout? (lit. OU réduction)
⎕← imprimer ce booléen en STDOUT
: et si oui:
  ⍴A forme d'une
   matrice de cette forme où chaque cellule est ses propres coordonnées
  =/¨ pour chaque cellule, les coordonnées verticales et horizontales sont-elles égales? (diagonale)
   multiplie les cellules de A par la
   matrice inverse (extrait la diagonale)
  D← dans D (pour D iagonal)
   inverser (revenir à la normale)
  A- soustraire de la
  ⍵+.×⍨ matrice A multiplier (la même chose que le produit scalaire, d'où le .) qu'avec x
  b- soustraire celle du
  D+.× produit de la matrice b de D et
  ⍺∇ appliquer cette fonction avec un A donné et que comme nouvelle valeur de x

Adam
la source
La sortie doit être le nombre d'itérations nécessaires pour une précision de e.
Zgarb
-1: Ce n'est pas une solution. Vous avez besoin de x0 car le but est de savoir combien d'étapes il faut pour atteindre la précision souhaitée à partir d'un point de départ particulier.
user1997744
@ user1997744 Ah, j'ai mal compris le problème. Pardon.
Adám
@ user1997744 Mieux?
Adám
1
@ user1997744 Pas d'opération arithmétique, juste la capacité de lire unaire , où en effet 0 n'est rien .
Adám
1

Python 3 , 132 octets

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

Essayez-le en ligne!

Utilise une solution récursive.

notjagan
la source
@ Adám, je ne suis pas sûr de bien comprendre. Je l'ai interprété comme fn'ayant pas le nom dans le bloc de code, que j'ai maintenant corrigé; cependant, s'il s'agit d'un problème entièrement différent, il peut toujours s'agir d'un problème.
notjagan
@ Adám Cette réponse semble corroborer ce que j'ai actuellement; c'est une fonction avec du code d'aide qui est capable de fonctionner comme une unité après sa définition.
notjagan
Ah ok. Pas de soucis alors. Je ne connais pas Python, donc j'étais simplement curieux. Bon travail!
Adám
Le critère d'arrêt n'est-il pas "Cela signifie que chaque composante des vecteurs en magnitude absolue diffère d'au plus e"? Fondamentalement, la norme max, pas la norme L2.
NikoNyrh
@NikoNyrh Fixed.
notjagan
1

R , 138 octets

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}

Essayez-le en ligne!

merci à NikoNyrh pour avoir corrigé un bug

Il convient également de noter qu'il existe un package R, Rlinsolvequi contient une lsolve.jacobifonction, renvoyant une liste avec x(la solution) et iter(le nombre d'itérations nécessaires), mais je ne suis pas sûr qu'il effectue les calculs corrects.

Giuseppe
la source
Le critère d'arrêt n'est-il pas "Cela signifie que chaque composante des vecteurs en magnitude absolue diffère d'au plus e"? Fondamentalement, la norme max, pas la norme L2.
NikoNyrh
@NikoNyrh vous avez raison! heureusement, la normfonction fournit cela pour moi aussi sans octets supplémentaires.
Giuseppe
1

Clojure, 212 198 196 octets

#(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)))

Implémenté sans bibliothèque matricielle, il itère le processus 1e9 fois pour obtenir la bonne réponse. Cela ne fonctionnerait pas sur des entrées trop mal conditionnées mais devrait fonctionner correctement dans la pratique.

Moins golfé, j'étais assez satisfait des expressions de Ret D:) La première entrée %(A) doit être un vecteur, pas une liste pour assocpouvoir être utilisée.

(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
la source