I’m having a problem trying to run my code. I am getting the message that the index is out of bounds. I tried running the part that has being showing problem separately for n = 9, c = 4 and U with dimision 1x81, as you can see bellow, it seems to be working
U = [-8.41470985e-01, -6.81638760e-01, -4.79425539e-01, -2.47403959e-01,
3.33066907e-16, 2.47403959e-01, 4.79425539e-01, 6.81638760e-01,
8.41470985e-01, -8.33843761e-01, -8.00032668e-01, -6.08192095e-01,
-3.25207319e-01, 9.80459808e-15, 3.25207319e-01, 6.08192095e-01,
8.00032668e-01, 8.33843761e-01, -5.81147313e-01, -5.76435639e-01,
-4.47982795e-01, -2.41978646e-01, 1.36309486e-14, 2.41978646e-01,
4.47982795e-01, 5.76435639e-01, 5.81147313e-01, -1.19943972e-01,
-1.16307504e-01, -8.53355534e-02,-4.34878176e-02, 9.57975558e-15,
4.34878176e-02, 8.53355534e-02, 1.16307504e-01, 1.19943972e-01,
3.78965423e-01, 3.89177337e-01, 3.19770624e-01, 1.80542815e-01,
7.10981512e-18, -1.80542815e-01, -3.19770624e-01,-3.89177337e-01,
-3.78965423e-01, 7.15338958e-01, 7.31044969e-01, 5.94841259e-01,
3.33049195e-01, -9.32965702e-15, -3.33049195e-01, -5.94841259e-01,
-7.31044969e-01, -7.15338958e-01, 7.52048363e-01, 7.67919215e-01,
6.23724461e-01, 3.48656958e-01, -1.30101869e-14, -3.48656958e-01,
-6.23724461e-01, -7.67919215e-01, -7.52048363e-01, 4.74997926e-01,
4.84909559e-01, 3.93652622e-01 , 2.19943579e-01, -9.08171401e-15,
-2.19943579e-01, -3.93652622e-01, -4.84909559e-01, -4.74997926e-01,
-8.88178420e-16, 4.44089210e-15, -2.55351296e-15, 3.33066907e-15,
-8.60422844e-16, 1.88737914e-15, -8.88178420e-16, -2.22044605e-16,
1.99840144e-15]
n = 9
N = 3
c = (n-1)/(N-1)
U_aux = []
for i in range(n):
for j in range(n):
k = i*n + j
if (i%c == 0) and (j%c == 0):
U_aux.append(U[k])
print(U_aux, '\n')
But when I came back to my code it isn’t. Since it shows the error: “index 29403 is out of bounds for axis 0 with size 243”. Can someone help me please?
import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.sparse import csc_matrix
from scipy.sparse import lil_matrix
import scipy.sparse.linalg as spla
def ij2n(i,j,n):
np_= i+j*n
return np_
a = -1
b = 1
def h_(n):
h = (b-a)/(n-1)
return h
def SolveSystem(n):
A = lil_matrix((n**2, n**2))
for i in range(n):
for j in range(n):
np_ = ij2n(i,j,n)
np_right = ij2n(i+1,j,n)
np_left = ij2n(i-1,j,n)
np_top = ij2n(i,j+1,n)
np_bottom = ij2n(i,j-1,n)
if j==0 or j == n-1:
A[np_,np_] = 1
elif i != 0 and i != n-1:
A[np_,np_] = (-4+10*pow(h,2))/pow(h,2)
A[np_,np_right] = 1/pow(h,2)
A[np_,np_left] = 1/pow(h,2)
A[np_,np_top] = 1/pow(h,2)
A[np_,np_bottom] = 1/pow(h,2)
elif i == 0:
A[np_,np_] = (-4+10*pow(h,2)-h)/pow(h,2)
A[np_,np_right] = 2/pow(h,2)
A[np_,np_top] = 1/pow(h,2)
A[np_,np_bottom] = 1/pow(h,2)
elif i == n-1:
A[np_,np_] = (-4+10*pow(h,2)-h)/pow(h,2)
A[np_,np_left] = 2/pow(h,2)
A[np_,np_top] = 1/pow(h,2)
A[np_,np_bottom] = 1/pow(h,2)
#print(A,'\n')
B1 = []
B2 = []
x = []
for i in range (pow(n,2)):
if i<=n-1:
x.append(-1+(i*h))
B1.append(np.sin(x[i]))
else:
B2.append(0)
B = np.transpose(np.concatenate((B1, B2)))
#Resolution of the system AU=B
A_csc = csc_matrix(A)
U = spla.spsolve(A_csc, B)
return U
# Setting number of point for the solution
n = 243
h = h_(n)
U = SolveSystem(n)
#Defining x_0, x_f, y_0 and y_f
x_0 = -1
x_f = 1
y_0 = -1
y_f = 1
# plotting
x = np.linspace(x_0,x_f,n)
y = np.linspace(y_0,y_f,n)
Xref, Yref = np.meshgrid(x,y)
U=np.reshape(U,Xref.shape)
fig = plt.figure(figsize = (8,8))
ax = plt.axes(projection='3d')
ax.plot_surface(Xref, Yref, U, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');
# rotation
ax.view_init(30, 60)
plt.savefig('3Dplot1.png',dpi=600)
plt.show()
N = [3, 9, 27, 81]
hvals = []
for i in range(4):
hvals.append(h_(N[i]))
U_aux0 = SolveSystem(N[0])
U_aux1 = SolveSystem(N[1])
U_aux2 = SolveSystem(N[2])
U_aux3 = SolveSystem(N[3])
def function(x):
c = (n-1)/(N[x]-1)
U_aux = []
for i in range(n):
for j in range(n):
k = i*n + j
if (i%c == 0) and (j%c == 0):
U_aux.append(U[k])
return U_aux
E = np.zeros(4)
U_aux = function(0)
E[0] = np.subtract(U_aux-U_aux0)
print(E[0],'\n')
U_aux = function(1)
E[1] = = np.subtract(U_aux-U_aux1)
U_aux = function(2)
E[2] = = np.subtract(U_aux-U_aux2)
U_aux = function(3)
E[3] = = np.subtract(U_aux-U_aux3)