Resolva uma equação de matriz pelo método de Jacobi (revisado)

11

Formação Matemática

Seja A uma matriz N de N de números reais, vetor ba de N números reais e xa vetor N números reais desconhecidos. Uma equação da matriz é Ax = b.

O método de Jacobi é o seguinte: decomponha A = D + R, onde D é a matriz de diagonais e R são as entradas restantes.

se você fizer uma solução inicial de adivinhação x0, uma solução aprimorada será x1 = inversa (D) * (b - Rx), em que todas as multiplicações são multiplicação de vetores matriciais e inversa (D) é inversa da matriz.


Especificação do Problema

  • Entrada : Seu programa completo deve aceitar como entrada os seguintes dados: a matriz A, o vetor b, um palpite inicial x0 e um número de 'erro' e.
  • Saída : O programa deve gerar o número mínimo de iterações, de modo que a solução mais recente seja diferente da solução verdadeira, no máximo e. Isso significa que cada componente dos vetores em magnitude absoluta difere no máximo e. Você deve usar o método de Jacobi para iterações.

Como os dados são inseridos é sua escolha ; pode ser sua própria sintaxe em uma linha de comando; você pode receber informações de um arquivo, o que quiser.

Como os dados são gerados é sua escolha ; pode ser gravado em um arquivo, exibido na linha de comando, gravado como arte ASCII, qualquer coisa, desde que seja legível por um humano.

Detalhes adicionais

Você não recebe a solução verdadeira: como você calcula a solução verdadeira depende inteiramente de você. Você pode resolvê-lo pela regra de Cramer, por exemplo, ou computando um inverso diretamente. O que importa é que você tenha uma solução verdadeira para poder comparar com as iterações.

Precisão é um problema; as 'soluções exatas' de algumas pessoas para comparação podem ser diferentes. Para os fins deste código de golfe, a solução exata deve ser verdadeira com 10 casas decimais.

Para ser absolutamente claro, se mesmo um componente da sua solução de iteração atual exceder o componente correspondente na solução verdadeira por e, será necessário continuar iterando.

O limite superior de N varia de acordo com o hardware que você está usando e quanto tempo está disposto a gastar executando o programa. Para os fins deste código de golfe, assuma o máximo de N = 50.

Condições prévias

Quando seu programa é chamado, você pode assumir que o seguinte é válido o tempo todo:

  • N> 1 e N <51, ou seja, você nunca receberá uma equação escalar, sempre uma equação matricial.
  • Todas as entradas estão no campo de números reais e nunca serão complexas.
  • A matriz A é sempre tal que o método converge para a solução verdadeira, de modo que você sempre pode encontrar várias iterações para minimizar o erro (conforme definido acima) abaixo ou igual a e.
  • A nunca é a matriz zero ou a matriz de identidade, ou seja, existe uma solução.

Casos de teste

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

A verdadeira solução é (0,586, 1,138). A primeira iteração fornece x1 = (5/9, 1), diferindo em mais de 0,04 da solução verdadeira, em pelo menos um componente. Fazendo outra iteração, encontramos x2 = (0,555, 1,148), que difere em menos de 0,04 de (0,586, 1,138). Assim, a saída é

2

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

Nesse caso, a verdadeira solução é (2.2, -0.8) e o palpite inicial x0 já possui um erro menor que e = 1.0, portanto, produzimos 0. Ou seja, sempre que você não precisar fazer uma iteração, basta gerar

0

Avaliação da submissão

Este é o código de golfe, com todas as brechas padrão aqui proibidas. O menor programa completo correto (ou função), ou seja, vence o menor número de bytes. É desencorajado o uso de coisas como o Mathematica, que agrupam muitas das etapas necessárias em uma função, mas usam o idioma desejado.

user1997744
fonte
2
Você realmente deve esperar para obter mais feedback, especialmente considerando o recente post fechado. Os desafios do PPCG geralmente compartilham estrutura comum nas especificações, o que geralmente contribui para que sejam fáceis de entender, em vez de cansativos e ambíguos. Tente analisar alguns dos desafios razoavelmente votados e imitar o formato.
Uriel
@ Uriel Entendo isso, mas acho que fui exaustivo em minhas especificações, e o formato, embora não corresponda exatamente à maioria das perguntas, pode ser lido linearmente e orientar o leitor de acordo. O formato também deve ter em mente o conteúdo do problema.
user1997744
3
"O programa completo correto mais curto " parece que você só permite programas e não funções. Eu adicionaria "/ function".
Adám 15/11/19
2
+1 marcas de formatação ou quebra a capacidade do meu cérebro para se concentrar em uma questão
Stephen
1
@ user1997744 Sim, faz sentido. Acredito que o padrão é que qualquer outro código, como outras funções ou importações de python, seja permitido, mas também incluído no bytecount.
Stephen

Respostas:

4

APL (Dyalog) , 78 68 65 49 bytes

Exatamente o tipo de problema para o qual o APL foi criado.

-3 graças a Erik, o Outgolfer . -11 graças a ngn .

Função de infixo anônimo. Toma A como argumento à esquerda ex como argumento à direita. Imprime o resultado em STDOUT como unário vertical usando 1como marcas de contagem, seguidas por 0como pontuação. Isso significa que mesmo um resultado 0 pode ser visto, não sendo 1s antes do 0.

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

Experimente online!

Explicação na ordem de leitura

Observe como o código lê de maneira muito semelhante à especificação do problema:

{... } nos dados A, bec, e no x fornecido,
⎕← imprima
∨/ se existe alguma verdade na afirmação de que
e< e é menor que
|⍵- o valor absoluto de x menos
b⌹ b matriz-dividido pelo
⊃A b e primeiro de A, b e e (ie A),
←⍺ que é o argumento à esquerda
: e, em caso afirmativo,
  ⍺∇ repetir em
  D+.× D vezes da matriz
  b- b menos
  ⍵+.×⍨ x, matriz multiplicada por
  A- A menos
  ⌹D o inverso de D (as entradas restantes) em
   que D é
   A onde
  =/¨ existem
   coordenadas iguais para
  ⍴A a forma de A (ou seja, a diagonal)

Explicação passo a passo

A ordem real de execução da direita para a esquerda:

{} Função anônima onde A é e ⍵ é x:
A b c←⍺ divida o argumento esquerdo em A, bec
 escolha a primeira (A)
b⌹ divisão matricial com b (fornece o valor verdadeiro de x)
⍵- diferenças entre os valores atuais de x e os
| valores absolutos
e< aceitáveis erro menor do que aqueles?
∨/ verdadeiro para qualquer? (lit. OR redução)
⎕← imprime esse booleano para STDOUT
: e, se sim:
  ⍴A forma de Uma
   matriz dessa forma em que cada célula tem suas próprias coordenadas
  =/¨ para cada célula, as coordenadas verticais e horizontais são iguais? (diagonal)
   multiplique as células de A pelo  armazenamento
   inverso da matriz that (extrai a diagonal)
  D←em D (para D iagonal)
  
  A- subtraia  inversa (de volta ao normal) de uma
  ⍵+.×⍨ matriz A multiplique (a mesma coisa que produto escalar, daí a .) que com x
  b- subtraia aquela de b
  D+.× produto de matriz de D e que
  ⍺∇ aplique essa função com A dado e seja como novo valor de x

Adão
fonte
A saída deve ser o número de iterações necessárias para uma precisão de e.
Zgarb
-1: Esta não é uma solução. Você precisa de x0, pois o ponto principal é saber quantas etapas são necessárias para obter a precisão desejada a partir de um ponto de partida específico.
user1997744
@ user1997744 Ah, eu entendi errado o problema. Desculpe.
Adám 15/11
@ user1997744 Melhor?
Adám 15/11
1
@ user1997744 Não é operação aritmética, apenas capacidade de ler unário , onde de fato 0 é nada .
Adám 15/11/19
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

Experimente online!

Usa uma solução recursiva.

notjagan
fonte
@ Adám Não tenho certeza se entendi bem. Eu interpretei como fnão tendo o nome dentro do bloco de código, que eu corrigi agora; no entanto, se for um problema completamente diferente, ainda poderá ser um problema.
notjagan
@ Adám Essa resposta parece corroborar o que tenho atualmente; é uma função com código auxiliar capaz de funcionar como uma unidade após sua definição.
notjagan
Ah ok. Deixa pra lá então. Eu não conheço Python, então fiquei curioso. Bom trabalho!
Adám
O critério de parada não é "Isso significa que cada componente dos vetores em magnitude absoluta difere no máximo e"? Basicamente, a norma máxima, não a norma L2.
NikoNyrh
@NikoNyrh Fixed.
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}

Experimente online!

graças ao NikoNyrh por corrigir um erro

Também vale a pena notar que existe um pacote R, Rlinsolve que contém uma lsolve.jacobifunção, retornando uma lista com x(a solução) e iter(o número de iterações necessárias), mas não tenho certeza de que ele faça os cálculos corretos.

Giuseppe
fonte
O critério de parada não é "Isso significa que cada componente dos vetores em magnitude absoluta difere no máximo e"? Basicamente, a norma máxima, não a norma L2.
NikoNyrh
@NikoNyrh você está correto! felizmente, a normfunção fornece isso para mim e também para nenhum bytes adicional.
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 sem uma biblioteca de matrizes, itera o processo 1e9 vezes para obter a resposta correta. Isso não funcionaria em entradas muito mal condicionadas, mas deveria funcionar bem na prática.

Menos golfista, fiquei bastante satisfeito com as expressões de Re D:) A primeira entrada %(A) deve ser um vetor, não uma lista para que assocpossa ser usada.

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