Matriz exponencial de uma matriz assimétrica real com Fortran 95 e LAPACK
Recentemente, fiz uma pergunta na mesma linha para matrizes enviesadas-hermitianas. Inspirado pelo sucesso dessa pergunta, e depois de bater minha cabeça contra uma parede por algumas horas, estou olhando para o exponencial da matriz de matrizes assimétricas reais. O caminho para encontrar os autovalores e autovetores parece bastante complicado, e receio ter me perdido.
Antecedentes: Há algum tempo, fiz essa pergunta sobre a física teórica SE. O resultado me permite formular equações mestres como matrizes assimétricas reais. No caso independente do tempo, a equação principal é resolvida exponenciando essa matriz. No caso dependente do tempo, será necessária integração. Só estou preocupado com a independência do tempo no momento.
Ao olhar para as várias sub-rotinas eu acho que deveria ser chamado ( ? Gehrd , ? Orghr , ? Hseqr ...) não está claro se ele seria mais simples para lançar a matriz de real*8a complex*16e prossiga com as versões duplas complexos destas rotinas, ou fique com real*8e tente dobrar o número de minhas matrizes e criar uma matriz complexa delas mais tarde.
Então, quais rotinas devo chamar (e em que ordem) e devo usar as versões duplas reais ou as versões duplas complexas? Abaixo está uma tentativa de fazer isso com versões duplas reais. Eu fiquei preso encontrando os autovalores e autovetores de L*t.
Gostaria de pensar seriamente se a matriz é realmente completamente arbitrária: existe alguma transformação que a tornaria hermitiana? A física garante que a matriz seja diagonalizável (com uma matriz autovetor razoavelmente condicionada)?
Se realmente não houver simetria a ser explorada, comece lendo Dezenove maneiras duvidosas de calcular o exponencial da matriz , que é a referência padrão (e é escrita pelo autor do MATLAB e co-autor do G & vL) .
O melhor que posso fazer é transformá-lo em uma matriz diagonal de blocos com blocos assimétricos. Isso por si só é muito interessante. A maioria desses blocos é , e eu posso resolvê-los analiticamente. Existe um bloco restante sem simetrias para explorar. 2×22×24×44×4
Eu gosto desta resposta; o caso assimétrico possui armadilhas suficientes que vale a pena considerar se pode haver uma formulação do seu problema que leve a matrizes simétricas em vez de matrizes assimétricas.
@ MarkS.Everitt: Você parece estar quase lá ... qual o tamanho das matrizes? ~ 36 x 36 de novo?
— 21412 Jack Poulson
Nesse caso, , mas existe a possibilidade de subir para . 16×1616×1636×3636×36
@ MarkS.Everitt: Então, o seu problema agora é exatamente como exponenciar matrizes 4x4. Isso é suficientemente pequeno para que a análise assintótica seja irrelevante; portanto, a resposta dependerá inteiramente dos valores. Eu realmente não posso dizer mais, a menos que você traduza sua publicação de física vinculada em álgebra linear (o que é um superoperador?!?).
Para desenvolver o que Jack disse, a abordagem padrão que parece ser usada em software (como EXPOKIT, mencionado em sua pergunta anterior) é escalar e esquadrar, seguida pela aproximação de Padé (métodos 2 e 3) ou pelos métodos do subespaço de Krylov (método 20) Em particular, se você estiver olhando para integradores exponenciais, considere os métodos do subespaço de Krylov e veja artigos sobre integradores exponenciais (algumas referências são mencionadas juntamente com o Método 20 no artigo de Moler e van Loan).
Se você está decidido a usar vetores próprios, considere o uso de sistemas triangulares de vetores próprios (método 15); como sua matriz pode ser não-diagonalizável, essa abordagem pode não ser a melhor opção, mas é melhor do que tentar calcular os vetores próprios e os valores próprios diretamente (ou seja, método 14).
Reduzir a forma de Hessenberg não é uma boa ideia (método 13).
Não está claro para mim se você seria melhor atendido com aritmética real ou complexa, uma vez que a aritmética complexa de Fortran é rápida, mas pode estourar / estourar (consulte "Qual é realmente melhor os compiladores Fortran?" ).
Você pode ignorar com segurança os Métodos 5-7 (os métodos baseados no solucionador de EDO são ineficientes), os Métodos 8 a 13 (caro), o Método 14 (o cálculo de vetores próprios de matrizes grandes é difícil, sem estrutura especial e propenso a erros numéricos em casos mal condicionados). e Método 16 (calcular a decomposição de Jordan de uma matriz é numericamente instável). Os métodos 17-19 são mais difíceis de implementar; em particular, os métodos 17 e 18 exigiriam mais leitura. O método 1 é uma opção de retorno para escalar e esquadrar se as aproximações de Padé não funcionarem bem.
Editar # 1 : Com base nos comentários em resposta à resposta de Jack, a diagonalização de blocos parece ser uma opção; nesse caso, algo como o Método 18 (diagonalização triangular de blocos) é um método muito bom de usar. Hesitei em recomendá-lo no início porque sua pergunta não especificou essa estrutura, mas se você tiver uma transformação que bloqueie a diagonalização de sua matriz, isso tira a maior parte da complexidade da abordagem. Você só quer ter certeza de usar o truque de GW Stewart para decompor cada matriz diagonal de bloco emBjBj
onde é a média dos valores próprios da matriz diagonal do ésimo bloco. Essa decomposição tornará quase nilpotente, o que aumentará a precisão dos seus cálculos exponenciais da matriz. Esse truque é discutido na página 26 da versão do artigo de Moler & van Loan ao qual Jack vinculou.γjγjjjEjEj
Votado, mas o pessoal que implementou o LAPACK não é ingênuo em relação à aritmética complexa, especialmente considerando o tempo todo que Kahan dedicou a analisá-la. No entanto, ele ainda deve computar na aritmética real se ela salvar flops; o custo da conversão posteriormente é apenas versus toda a álgebra linear densa de . O(n2)O(n2)O(n3)O(n3)
Sem dúvida, eles sabem o que estão fazendo; Não estou preocupado com a implementação do LAPACK. Estou mais surpreso com o comportamento do compilador Fortran.
Sim, o compilador pode ser mais um problema do que o LAPACK bem escrito. Ele pode ser desconcertante para descobrir que seu programa só foi falhando porque as implementações de valor absoluto e divisão utilizada pelo compilador foram mal feita ...
Eu tenho uma sub-rotina Fortran simples que calcula o expoente de uma matriz arbitrária. Eu verifiquei o comando Matlab e está tudo bem. Baseia-se em escala e quadratura. Eu escrevi alguns anos atrás.
Eu gostaria de definir outra sub-rotina, como as que eu baixo de gams.nist.gov. Mas sem sorte ainda.