Skip to main content
1 de 2

Como criar uma matriz tridiagonal em Python?

Olá, eu gostaria de saber como escrever uma matriz tridiagonal NxN. Eu tenho três vetores que serão atualizados em um loop.

Eu estou trabalhando com algo semlhnante a isso: Example

Nota: Eu apenas gostaria de saber sobre zeros e que parâmetros do Python eu poderia usar para ajustar isso. Bem, eu tenho dois códigos aqui, o primeiro eu escrevi em Fortran e está funcionando bem. O segundo é o que eu tentei escrever em Python.

Fortran:

do i=2,n-1
  do j=2,n-1
    if (i.eq.j) then

      D(i,j+1)=-u_med(i+1)/(delta_r(i)*delta_r(i+1))

      t1 =u_med(i+1)/(delta_r(i)*delta_r(i))

      t2 = u_med(i)/(delta_r(i)*delta_r(i))

      D(i,j)= t1 + t2 + V(i)

      D(i+1,j)=-u_med(i+1)/(delta_r(i)*delta_r(i+1))

    end if

  end do
end do

Python:

for i in range(2,n):
  for j in range(2,n):
    if i == j:
      D[i][j+1] = - u_med[i+1]/(delta_r[i]*delta_r[i+1])
      t1 = u_med[i+1]/(delta_r[i]*delta_r[i])
      t2 = u_med[i]/(delta_r[i]*delta_r[i])
      D[i][j]= t1 + t2 + V[i]
      D[i+1][j]= - u_med[i+1]/(delta_r[i]*delta_r[i+1])

t1 = u_med[2]/(delta_r[1]*delta_r[1])
t2 = 0
D[1][1]= t1 + t2 + V[1]
D[1][2]= - u_med[2]/(delta_r[1]*delta_r[2])
D[2,1]= - u_med[2]/(delta_r[2]*delta_r[1])
t1 = 0
t2 = u_med[n]/(delta_r[n]*delta_r[n])

D[n][n]= t1 + t2 + V[n]

Que dá o erro:

D[i][j+1] = - u_med[i+1]/(delta_r[i]*delta_r[i+1]) ValueError: setting an array element with a sequence.

Definições de u_med, delta_r, r, n e D em Python:

u_med = np.zeros((n,2))
for i in range(2,n):
  tta1 = r[i]*u[i]
  tta2 = r[i-1]*u[i-1]
  u_med[i] = 0.5*(tta1 + tta2)/(r[i] - r[i-1])

u_med[1] = u_med[2]

delta_r = np.zeros((n-1,2))
for i in range(2,n-1):
  ft1 = r[i+1]*r[i+1]
  ft2 = r[i-1]*r[i-1]
  ft3 = 2*r[i]*(r[i+1] - r[i-1])
  delta_r[i] = math.sqrt(0.125*abs(ft1 - ft2 + ft3))

ri=0
n1 = 51

r1 = ri
r2 = 250
hr1 = (r2-r1)/(n1-1)

r = np.zeros((n1,1))
for i in range(n1):
  r[i] = r1 + i*hr1

u = np.zeros((n+1,1))
for i in range(1,n+1):
  i = 1

D = npm.zeros((n,n))

Eu suspeito que o erro esteja em:

D = npm.zeros((n,n))

Mas como eu poderia reescrever isso?