Resolver una ecuación matricial por el método de Jacobi (Revisado)

11

Antecedentes matemáticos

Sea A una matriz N de N de números reales, un vector ba de N números reales y un vector N de números reales desconocidos. Una ecuación matricial es Ax = b.

El método de Jacobi es el siguiente: descomponer A = D + R, donde D es la matriz de diagonales y R son las entradas restantes.

Si realiza una solución inicial de conjetura x0, una solución mejorada es x1 = inversa (D) * (b - Rx) donde todas las multiplicaciones son multiplicación de matriz-vector e inversa (D) es la matriz inversa.


Especificación del problema

  • Entrada : Su programa completo debe aceptar como entrada los siguientes datos: la matriz A, el vector b, una suposición inicial x0 y un número de 'error' e.
  • Salida : El programa debe generar la cantidad mínima de iteraciones de manera que la última solución difiera de la solución verdadera, a lo sumo e. Esto significa que cada componente de los vectores en magnitud absoluta difiere como máximo en e. Debe usar el método de Jacobi para las iteraciones.

Cómo se ingresan los datos es su elección ; podría ser su propia sintaxis en una línea de comando, podría tomar la entrada de un archivo, lo que elija.

Cómo se emiten los datos es su elección ; puede escribirse en un archivo, mostrarse en la línea de comando, escribirse como arte ASCII, cualquier cosa, siempre que sea legible por un humano.

Más detalles

No se le da la solución verdadera: la forma de calcular la solución verdadera depende completamente de usted. Puede resolverlo mediante la regla de Cramer, por ejemplo, o calculando un inverso directamente. Lo importante es que tenga una verdadera solución para poder comparar con las iteraciones.

La precisión es un problema; Las 'soluciones exactas' de algunas personas para la comparación pueden diferir. Para los propósitos de este código de golf, la solución exacta debe ser verdadera a 10 decimales.

Para ser absolutamente claro, si incluso un componente de su solución de iteración actual excede su componente correspondiente en la solución verdadera por e, entonces debe seguir iterando.

El límite superior a N varía según el hardware que esté utilizando y cuánto tiempo esté dispuesto a pasar ejecutando el programa. Para los propósitos de este código de golf, suponga un máximo de N = 50.

Condiciones previas

Cuando se llama a su programa, puede asumir que lo siguiente se cumple en todo momento:

  • N> 1 y N <51, es decir, nunca se le dará una ecuación escalar, siempre una ecuación matricial.
  • Todas las entradas están sobre el campo de números reales, y nunca serán complejas.
  • La matriz A siempre es tal que el método converge a la solución verdadera, de modo que siempre puede encontrar una serie de iteraciones para minimizar el error (como se definió anteriormente) a continuación o igual a e.
  • A nunca es la matriz cero o la matriz de identidad, es decir, hay una solución.

Casos de prueba

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

La verdadera solución es (0.586, 1.138). La primera iteración da x1 = (5/9, 1), que difiere en más de 0.04 de la solución verdadera, en al menos un componente. Tomando otra iteración encontramos, x2 = (0.555, 1.148) que difiere en menos de 0.04 de (0.586, 1.138). Por lo tanto, la salida es

2

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

En este caso, la verdadera solución es (2.2, -0.8) y la suposición inicial x0 ya tiene un error menor que e = 1.0, por lo tanto, arrojamos 0. Es decir, cuando no necesita hacer una iteración, simplemente genera

0

Evaluación de envío

Este es el código de golf, con todas las lagunas estándar por la presente no permitidas. El programa completo correcto más corto (o función), es decir, gana el menor número de bytes. Se desaconseja usar cosas como Mathematica que envuelven muchos de los pasos necesarios en una función, pero use el lenguaje que desee.

usuario1997744
fuente
2
Realmente deberías esperar para recibir más comentarios al respecto, especialmente dada la reciente publicación cerrada. Los desafíos de PPCG generalmente comparten una estructura común en las especificaciones, lo que generalmente contribuye a que sean fáciles de entender, en lugar de cansados ​​y ambiguos. Trate de ver algunos de los desafíos razonablemente votados e imite el formato.
Uriel
@Uriel Me doy cuenta de esto, pero creo que he sido exhaustivo en mi especificación, y el formato, aunque no se ajusta exactamente a la mayoría de las preguntas, puede leerse linealmente y guiar al lector en consecuencia. El formato también debe tener en cuenta el contenido del problema en sí.
user1997744
3
"El programa completo correcto más corto " suena como si solo permitiera programas y no funciones. Yo agregaría "/ function".
Adám
2
El formateo de +1 hace o rompe la capacidad de mi cerebro para concentrarse en una pregunta
Stephen
1
@ user1997744 Sí, tiene sentido. Creo que el valor predeterminado es que cualquier otro código, como otras funciones o importaciones de Python, está permitido, pero también incluido en el bytecount.
Stephen

Respuestas:

4

APL (Dyalog) , 78 68 65 49 bytes

Exactamente el tipo de problema para el que se creó APL.

-3 gracias a Erik the Outgolfer . -11 gracias a ngn .

Función de infijo anónimo. Toma A ser como argumento izquierdo y x como argumento derecho. Las impresiones dan como resultado STDOUT como unario vertical utilizando 1marcas de conteo, seguidas de 0puntuación. Esto significa que incluso se puede ver un resultado 0, no hay 1s antes del 0.

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

Pruébalo en línea!

Explicación en orden de lectura

Observe cómo el código se lee de manera muy similar a la especificación del problema:

{... } en la A, bye dada, yy la x dada,
⎕← imprima
∨/ si hay alguna verdad en la afirmación de que
e< e es menor que
|⍵- el valor absoluto de x menos
b⌹ b, dividido en matriz por
⊃A b e el primero de A, by e (es decir, A),
←⍺ que son el argumento izquierdo
: y, en caso afirmativo,
  ⍺∇ recurren en
  D+.× D matriz-veces
  b- b menos
  ⍵+.×⍨ x, matriz multiplicada por
  A- A menos
  ⌹D la inversa de D (las entradas restantes)
   donde D es
   A donde
  =/¨ hay
   coordenadas iguales para
  ⍴A la forma de A (es decir, la diagonal)

Explicación paso a paso

El orden real de ejecución de derecha a izquierda:

{} Función anónima donde es A be y ⍵ es x:
A b c←⍺ divide el argumento izquierdo en A, b, y e
 elige la primera
b⌹ división de matriz (A) con
⍵- diferencias b (da el valor verdadero de x) entre los valores actuales de x y esos
| valores absolutos
e< aceptables error menos que esos?
∨/ cierto para cualquiera? (encendido O reducción)
⎕← imprima ese booleano en STDOUT
: y, de ser así:
  ⍴A forma de una
   matriz de esa forma donde cada celda tiene sus propias coordenadas
  =/¨ para cada celda, ¿son iguales las coordenadas vertical y horizontal? (diagonal)
   multiplique las celdas de A con el  almacenamiento
   inverso de la matriz que (extrae diagonal)
  D←en D (para D iagonal)
   inversa (volver a la normalidad)
  A- restar de la
  ⍵+.×⍨ matriz A multiplicar (lo mismo que el producto de punto, de ahí el .) que con x
  b- reste eso del
  D+.× producto de la matriz b de D y que
  ⍺∇ aplique esta función con A dado y que como nuevo valor de x

Adán
fuente
La salida debe ser el número de iteraciones requeridas para una precisión de e.
Zgarb
-1: Esto no es una solución. Necesita x0 ya que el punto es saber cuántos pasos se necesitan para llegar a la precisión deseada desde un punto de partida particular.
user1997744
@ user1997744 Ah, entendí mal el problema. Lo siento.
Adám
@ user1997744 ¿Mejor?
Adám
1
@ user1997744 No es una operación aritmética, solo la capacidad de leer unario , donde de hecho 0 no es nada .
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

Pruébalo en línea!

Utiliza una solución recursiva.

notjagan
fuente
@ Adám, no estoy seguro de entenderlo. Lo interpreté como que fno tenía el nombre dentro del bloque de código, que ahora he arreglado; sin embargo, si se trata de un problema completamente diferente, aún puede ser un problema.
notjagan
@ Adám Esa respuesta parece corroborar lo que tengo actualmente; Es una función con código auxiliar que es capaz de funcionar como una unidad después de su definición.
notjagan
Ah ok Olvidalo entonces. No conozco Python, así que tenía curiosidad. ¡Buen trabajo!
Adám
¿No es el criterio de detención "Esto significa que cada componente de los vectores en magnitud absoluta difiere como máximo en e"? Básicamente la norma máxima, no la norma L2.
NikoNyrh
@NikoNyrh fijo.
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}

Pruébalo en línea!

gracias a NikoNyrh por arreglar un error

También vale la pena señalar que hay un paquete R, Rlinsolveque contiene una lsolve.jacobifunción, que devuelve una lista con x(la solución) y iter(el número de iteraciones requeridas), pero no estoy seguro de que haga los cálculos correctos.

Giuseppe
fuente
¿No es el criterio de detención "Esto significa que cada componente de los vectores en magnitud absoluta difiere como máximo en e"? Básicamente la norma máxima, no la norma L2.
NikoNyrh
@NikoNyrh tienes razón! Afortunadamente, la normfunción también proporciona eso para mí sin bytes adicionales.
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)))

Implementado sin una biblioteca matricial, itera el proceso 1e9 veces para obtener la respuesta correcta. Esto no funcionaría con entradas demasiado mal acondicionadas, pero debería funcionar bien en la práctica.

Menos golf, estaba bastante contento con las expresiones de Ry D:) La primera entrada %(A) tiene que ser un vector, no una lista para que assocpueda usarse.

(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
fuente