PDE Solving: NoneType не имеет атрибута «toarray», ошибка

Question

Я учусь использовать методы конечных разностей в Python, используя в качестве примера уравнение Пуассона. Но когда я запускаю все нормально, до строки кода A=ASparse.toarray() в обратном разделе, где я получаю ошибку, упомянутую в заголовке поста. Почему это происходит, откуда и как мне это решить?

Код:

import numpy as np
import scipy
import scipy.sparse as sps
import matplotlib.pyplot as plt


# Elliptic PDE problems can be converted to Linear Algebra problems!!!

# Begin with making the matrix A for conversion to a linear system

J = 5
L = 4
I = (J - 1) * (L - 1)

LeadDiag = -4 * np.ones(I)
diag=[LeadDiag]
ASparse = sps.diags(diagonals=diag, offsets=[0])
ASparse = sps.csc_matrix(ASparse)
A = ASparse.toarray()


UpperDiag = np.ones(I - 1)
for i in range(0, I - 1):
    if np.mod(i+1, L-1) == 0:
        UpperDiag[i] = 0

LowerDiag = np.ones(I - 1)
for i in range(0, I - 1):
    if np.mod(i+1, L-1)==0:
        LowerDiag[i]=0

SuperDiag = np.ones(I - (L - 1))
SubDiag = np.ones(I - (L - 1))

Diag = [LeadDiag, UpperDiag, LowerDiag, SuperDiag, SubDiag]
ASparse = sps.diags(diagonals=Diag, offsets=[0, 1, -1, L-1, -(L-1)])
ASparse = sps.csc_matrix(ASparse)
A = ASparse.toarray()
print(A)

# Short way to make A

def createA(J, L):
    I = (J - 1)*(L - 1)

LeadDiag = -4 * np.ones(I)

UpperDiag = np.ones(I - 1)
for i in range(0, I - 1):
    if np.mod(i + 1, L - 1)==0:
        UpperDiag[i] = 0

LowerDiag = np.ones(I - 1)
for i in range(1, I - 1):
    if np.mod(i+1, L-1)==0:
        LowerDiag[i] = 0

SuperDiag=np.ones(I - (L - 1))
SubDiag=np.ones(I-(L-1))

Diag=[LeadDiag, UpperDiag, LowerDiag, SuperDiag, SubDiag]
ASparse = sps.diags(diagonals=Diag, offsets=[0, 1, -1, L-1, -(L-1)])
ASparse=sps.csc_matrix(ASparse)

print(ASparse)

# Define f(x,y) (rho) and plot it on the grid given by (x_j, y_j), choosing J = 100, L = 100, L_x = L_y = 100

def funct(x, y, L_x, L_y):
    A = -1/ (np.square(2 * np.pi / L_x) + np.square(np.pi / L_y))
    return np.sin(2 * np.pi / L_x * x) * np.sin(np.pi/ L_y * y)

J = 100
L = 100
L_x = 10
L_y = 10

x = np.linspace(0, L_x, J+1, endpoint=True)
y = np.linspace(0, L_y, L+1, endpoint=True)
xv, yv = np.meshgrid(x, y)
print(x.shape, y.shape)

f = funct(xv, yv, L_x, L_y)
print(f.shape)

plt.imshow(f)
plt.colorbar()
plt.show()

# BOUNDARIES

print(f.shape)

Delt = L_x/J
fTemp = np.copy(f)
fTemp[1, :] = fTemp[1, :]+0
fTemp[L, :] = fTemp[L, :]+0
fTemp[:, 1] = fTemp[:, 1]+0
fTemp[:, J] = fTemp[:, J]+0

f1D = np.square(Delt)*fTemp[1:J, 1:L].reshape((J-1)*(L-1))
print(f1D.shape)

#Define the matrix A and compute its inverse

import scipy.sparse.linalg

ASparse=createA(J, L)

A=ASparse.toarray()
print(A.shape)

ASparseInv = sps.linalg.inv(ASparse)

# Solve the Linear Algebra problem!

u=np.dot(ASparseInv.toarray(), f1D)
print(u.shape)

u2D = u.reshape((J - 1), (L-1))
print(u2D.shape)

# Analytical plot

def uexactfunct(x, y, L_x, L_y):
    A=-1/np.square(2 * np.pi/L_x) + np.square(np.pi/L_y)
    return A*np.sin(2 * np.pi/L_x * x)*np.sin(np.pi/L_y * y)

# Analytical vs numerical plots
plt.imshow(u2D)
plt.colorbar()
plt.show()

plt.imshow(uexactfunct(xv, yv, L_x, L_y))
plt.colorbar()
plt.show()

plt.plot(x[1:-1], u2D[J/2+1, :], 'x', label='numerical')
plt.plot(x[1:-1], uexactfunct(x[1:-1], y[J/2+2], L_x, L_y), '-', label='analytical')
plt.legend()
plt.show()

plt.plot(y[1:-1], u2D[:, L/2 + 1], 'x', label='numerical')
plt.plot(y[1:-1], uexactfunct(x[L/2 + 2], y[1:-1], L_x, L_y), '-', label='analytical')
plt.legend()
plt.show()

Answer 1

Ваша createA() - функция ничего не возвращает. Таким образом, ASparse равен None в тот момент, который, очевидно, не имеет toarray() метода. Таким образом, вы захотите реализовать createA для возврата I , но даже это просто возвращает целое число.

Что вы действительно хотите сделать:

def createA(J, L):
    I = (J - 1) * (L - 1)
    LeadDiag = -4 * np.ones(I)
    diag = [LeadDiag]
    ASparse = sps.diags(diagonals=diag, offsets=[0])
    ASparse = sps.csc_matrix(ASparse)
    return ASparse

В графическом материале в конце сценария вы должны приводить число в целые числа.

plt.plot(x[1:-1], u2D[J/2+1, :], 'x', label='numerical')

должно быть

plt.plot(x[1:-1], u2D[int(J/2+1), :], 'x', label='numerical')

Я бы порекомендовал вам использовать отладчик, так что вы понимаете, что происходит в ваших сценариях. Особенно, когда вы работаете с матрицами и индексами. Множество мелочей, чтобы ошибиться:)