2011-04-30 23 views
11

Me gustaría crear una matriz tridiagonal de bloque a partir de tres numpy.ndarray. ¿Hay alguna forma (directa) de hacer eso en python?Bloque tridiagonal matrix python

¡Gracias de antemano!

Saludos

+1

¿Quieres que el resultado sea otra ndarray, o está abierta a la utilización de una matriz dispersa por el resultado? – talonmies

Respuesta

7

También puede hacer esto con matrices numpy "normales" a través de la indexación de lujo:

import numpy as np 
data = np.zeros((10,10)) 
data[np.arange(5), np.arange(5)+2] = [5, 6, 7, 8, 9] 
data[np.arange(3)+4, np.arange(3)] = [1, 2, 3] 
print data 

(podría reemplazar esas llamadas a np.arange con np.r_ si quería ser más concisa Por ejemplo, en lugar de., utilice data[np.r_[:3]+4, np.r_[:3]])

Esto produce:

[[0 0 5 0 0 0 0 0 0 0] 
[0 0 0 6 0 0 0 0 0 0] 
[0 0 0 0 7 0 0 0 0 0] 
[0 0 0 0 0 8 0 0 0 0] 
[1 0 0 0 0 0 9 0 0 0] 
[0 2 0 0 0 0 0 0 0 0] 
[0 0 3 0 0 0 0 0 0 0] 
[0 0 0 0 0 0 0 0 0 0] 
[0 0 0 0 0 0 0 0 0 0] 
[0 0 0 0 0 0 0 0 0 0]] 

Sin embargo, si usted va a utilizar matrices dispersas de todos modos, echar un vistazo a scipy.sparse.spdiags. (Tenga en cuenta que tendrá que anteponer datos falsos en sus valores de fila Si va a colocar los datos en una posición diagonal con un valor positivo (por ejemplo, los de 3 en la posición 4 en el ejemplo))

Como rápida ejemplo:

import numpy as np 
import scipy as sp 
import scipy.sparse 

diag_rows = np.array([[1, 1, 1, 1, 1, 1, 1], 
         [2, 2, 2, 2, 2, 2, 2], 
         [0, 0, 0, 0, 3, 3, 3]]) 
positions = [-3, 0, 4] 
print sp.sparse.spdiags(diag_rows, positions, 10, 10).todense() 

Este rendimientos:

[[2 0 0 0 3 0 0 0 0 0] 
[0 2 0 0 0 3 0 0 0 0] 
[0 0 2 0 0 0 3 0 0 0] 
[1 0 0 2 0 0 0 0 0 0] 
[0 1 0 0 2 0 0 0 0 0] 
[0 0 1 0 0 2 0 0 0 0] 
[0 0 0 1 0 0 2 0 0 0] 
[0 0 0 0 1 0 0 0 0 0] 
[0 0 0 0 0 1 0 0 0 0] 
[0 0 0 0 0 0 1 0 0 0]] 
+0

¡Gracias chicos! –

19

Con arrays numpy "regular", utilizando numpy.diag:

def tridiag(a, b, c, k1=-1, k2=0, k3=1): 
    return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3) 

a = [1, 1]; b = [2, 2, 2]; c = [3, 3] 
A = tridiag(a, b, c) 
0

Mi respuesta se basa en la respuesta de @ TheCorwoodRep. Acabo de publicarlo porque hice algunos cambios para hacerlo más modular, para que funcione para diferentes órdenes de matrices y también para cambiar los valores de k1, k2, , es decir, que decida dónde aparece la diagonal, se ocupará de la desbordamiento automático Al llamar a la función, puede especificar qué valores deben aparecer en las diagonales.

import numpy as np 
def tridiag(T,x,y,z,k1=-1, k2=0, k3=1): 
    a = [x]*(T-abs(k1)); b = [y]*(T-abs(k2)); c = [z]*(T-abs(k3)) 
    return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3) 

D=tridiag(10,-1,2,-1) 
2

respuesta @TheCorwoodRep 's en realidad se puede hacer en una sola línea. No hay necesidad de una función separada.

np.eye(3,3,k=-1) + np.eye(3,3)*2 + np.eye(3,3,k=1)*3 

Esto produce:

array([[ 2., 3., 0.], 
     [ 1., 2., 3.], 
     [ 0., 1., 2.]]) 
2

utilizar la función scipy.sparse.diags.

Ejemplo:

from scipy.sparse import diags 
import numpy as np 
# 
n = 10 
k = np.array([np.ones(n-1),-2*np.ones(n),np.ones(n-1)]) 
offset = [-1,0,1] 
A = diags(k,offset).toarray() 

Esto devuelve:

array([[-2., 1., 0., 0., 0.], 
     [ 1., -2., 1., 0., 0.], 
     [ 0., 1., -2., 1., 0.], 
     [ 0., 0., 1., -2., 1.], 
     [ 0., 0., 0., 1., -2.]])