Tridiagonal Matrix using Python

Question

I looked on the past questions on tridiagonals but none seem to be experiencing the problem i'm having. I'm trying to form a tridiagonal stiffness matrix for the non uniform Poisson equation using scipy.sparse.spdiags but do not seem to be receiving a matrix as a result.

def Poisson_Stiffness(x0):

N = len(x0) - 1 #THE AMOUNT OF ELEMENTS (NOT THE AMOUNT OF POINTS) x0, x1, ... , x_N

h = np.zeros(N)
a = np.zeros(N+1)
b = np.zeros(N)

for i in range(N):
    h[i] = x0[i+1] - x0[i] #Length of each nonuniform element

a[0] = 1/h[0]
for i in range(1,N):
    a[i] = 1/h[i] + 1/h[i-1] #Main Diagonal of stiffness matrix
a[N] = 1/h[N-1]

for i in range(N):
    b[i] = -1/h[i] #Upper and lower diagonal of stiffness matrix.


Tridiagonal_Data = np.array([[a],[b],[b]])
Positions = [0, 1, -1]

Stiffness_Matrix = scipy.sparse.spdiags(Tridiagonal_Data,Positions,N+1,N+1)

print Stiffness_Matrix

As a result from this, with x0 = [0,0.3,0.4,0.7,1]; I am getting the stiffness matrix as :

Jamess-MBP:Poisson jamesmalone$ python Poisson1d.py
(0, 0)  [  3.33333333  13.33333333  13.33333333   6.66666667   3.33333333]
(1, 0)  [-3.33333333 -10.         -3.33333333 -3.33333333]

My question is why does it come out like this and not in matrix form?

I have tried changing the data types to see if that was the problem but i would receive errors such as (for .toarray()):

Jamess-MBP:Poisson jamesmalone$ python Poisson1d.py
Traceback (most recent call last):
  File "Poisson1d.py", line 33, in <module>
    Poisson_Stiffness([0,0.3,0.4,0.7,1])
  File "Poisson1d.py", line 29, in Poisson_Stiffness
    Stiffness_Matrix = scipy.sparse.spdiags(Tridiagonal_Data,Positions,N+1,N+1).toarray()
  File "/Users/jamesmalone/anaconda/lib/python2.7/site-packages/scipy/sparse/base.py", line 637, in toarray
return self.tocoo().toarray(order=order, out=out)
  File "/Users/jamesmalone/anaconda/lib/python2.7/site-packages/scipy/sparse/coo.py", line 275, in toarray
B.ravel('A'), fortran)

RuntimeError: internal error: failed to resolve data types

Thanks in advance.

wflynny · Accepted Answer · 2015-06-25 19:50:23Z

2

Try using scipy.sparse.diags. I also cleaned up your code because you are not taking advantage of numpy's strengths (broadcasting) with those for loops. Also cleaned up some formatting according to PEP8:

from scipy.sparse import diags

x0 = np.array(x0)
N = len(x0) - 1

h = x0[1:] - x0[:-1]

a = np.zeros(N+1)
a[0] = 1/h[0]
a[1:-1] = 1/h[1:] + 1/h[:-1]
a[-1] = 1/h[-1]

b = -1/h

data = [a.tolist(), b.tolist(), b.tolist()]
positions = [0, 1, -1]

stiffness_matrix = diags(data, positions, (N+1, N+1))

print stiffness_matrix.toarray()

With x0 = [0, 0.3, 0.4, 0.7, 1], this yields

[[  3.33333333  -3.33333333   0.           0.           0.        ]
 [ -3.33333333  13.33333333 -10.           0.           0.        ]
 [  0.         -10.          13.33333333  -3.33333333   0.        ]
 [  0.           0.          -3.33333333   6.66666667  -3.33333333]
 [  0.           0.           0.          -3.33333333   3.33333333]]

answered Jun 25, 2015 at 19:50

wflynny

18.6k6 gold badges50 silver badges69 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

James Over a year ago

Thanks alot, a great help :) I'm also a big fan of the broadcasting, wasn't aware i could do it that way before.

Collectives™ on Stack Overflow

Tridiagonal Matrix using Python

1 Answer 1

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related