Calculating a matrix of 144x144 using multiple loops

I have a problem calculating the kuu_s matrix ,which is a size of 144x144; the program is only calculating 121x121 elements leaving the first row and columns to 0, and the last 22x22 elements also resulting in 0.

import numpy as np
import scipy.integrate as spi
# Define the integrand function
def f0(x,n):
if n == 1:
return (1-x)
elif n == 2:
return x
else:
return np.sin((n - 2) * np.pi * x)
def f1(x,n):
if n == 1:
return -1
elif n == 2:
return 1
else:
return np.cos((n - 2)*np.pi*x)*(n-2)*np.pi
def f00(x, m, n):
return f0(x,m)*f0(x,n)
def f11(x, m, n):
return f1(x,m)*f1(x,n)
def f10(x, m, n):
return f1(x,m)*f0(x,n)
def f01(x, m, n):
return f0(x,m)*f1(x,n)
# Define the dimensions of the tensor
NMHT_0 = 12
# Create the tensor
I00 = I11 = I10 = I01 = np.zeros((NMHT_0, NMHT_0))
# Perform the integration
for m in range(1, NMHT_0 ):
for n in range(1, NMHT_0 ):
result1, error = spi.quad_vec(f00, 0, 1, args=(m, n)) result2, error2 = spi.quad_vec(f11, 0, 1, args=(m, n)) result3, error3 = spi.quad_vec(f10, 0, 1, args=(m, n)) result4, error4 = spi.quad_vec(f01, 0, 1, args=(m, n))
I00[m, n] = result1 I11[m, n] = result2  I10[m, n] = result3 I01[m, n] = result4
a = 1
b = 1
A = np.random.rand(5,5)
Kuu_s = np.zeros((NMHT_0 ** 2, NMHT_0 ** 2))
k = 0
for i in range(1, NMHT_0):
for j in range(1, NMHT_0):
q = 0
k += 1
for m in range(1, NMHT_0):
for n in range(1, NMHT_0):
q += 1
Kuu_s[k, q] += (
(A[1, 1] * (1 / a ** 2) * (I11[i, m] * I00[j, n])) +
(A[1, 3] * (1 / (a * b)) * (I10[i, m] * I01[j, n] + I01[i, m] * I10[j, n])) +
(A[3, 3] * (1 / b ** 2) * (I00[j, n] * I11[i, m]))
)
np.savetxt("matrixkuu2.docx",Kuu_s,delimiter='  ',fmt='%1.3e')

the following pic is from mathcad

I've tried changing the limits of m n i j ,but I get the error of out-of-bound index.
How could I fix the issue so that the program calculates all the 144x144 elements and adds kuu_s to each loop

I've fixed the spacing problems with your code, plus the issues raised in the comments. This is the result. It does produce a 144x144 output array.

import numpy as np
import scipy.integrate as spi
# Define the integrand function
def f0(x,n):
if n == 1:
return (1-x)
elif n == 2:
return x
else:
return np.sin((n - 2) * np.pi * x)
def f1(x,n):
if n == 1:
return -1
elif n == 2:
return 1
else:
return np.cos((n - 2)*np.pi*x)*(n-2)*np.pi
def f00(x, m, n):
return f0(x,m)*f0(x,n)
def f11(x, m, n):
return f1(x,m)*f1(x,n)
def f10(x, m, n):
return f1(x,m)*f0(x,n)
def f01(x, m, n):
return f0(x,m)*f1(x,n)
# Define the dimensions of the tensor
NMHT_0 = 12
# Create the tensor
I00 = np.zeros((NMHT_0, NMHT_0))
I11 = np.zeros((NMHT_0, NMHT_0))
I10 = np.zeros((NMHT_0, NMHT_0))
I01 = np.zeros((NMHT_0, NMHT_0))
# Perform the integration
for m in range(NMHT_0):
for n in range(NMHT_0):
result1, error1 = spi.quad_vec(f00, 0, 1, args=(m, n))
result2, error2 = spi.quad_vec(f11, 0, 1, args=(m, n))
result3, error3 = spi.quad_vec(f10, 0, 1, args=(m, n))
result4, error4 = spi.quad_vec(f01, 0, 1, args=(m, n))
I00[m, n] = result1 
I11[m, n] = result2  
I10[m, n] = result3 
I01[m, n] = result4
a = 1
b = 1
A = np.random.rand(5,5)
Kuu_s = np.zeros((NMHT_0 ** 2, NMHT_0 ** 2))
for i in range(NMHT_0):
for j in range(NMHT_0):
k = i * NMHT_0 + j
for m in range(NMHT_0):
for n in range(NMHT_0):
q = m * NMHT_0 + n
Kuu_s[k, q] += (
(A[1, 1] * (1 / a ** 2) * (I11[i, m] * I00[j, n])) +
(A[1, 3] * (1 / (a * b)) * (I10[i, m] * I01[j, n] + I01[i, m] * I10[j, n])) +
(A[3, 3] * (1 / b ** 2) * (I00[j, n] * I11[i, m]))
)
np.savetxt("matrixkuu2.txt",Kuu_s,delimiter='  ',fmt='%1.3e')

When creating Kuu_s, you start your loops at 1, when they should start at 0 because Python is zero-indexed. The second issue is starting your k and q values. Starting them at 0 and incrementing at the start of the loop means they will first be used when they equal 1, but (as I already mentioned) Pyhton is zero-indexed, so you're missing the first row and column. You have two options:

1. Start k and q at 0 and increment them at the end of the loop.
2. Start k and q at -1 and increment them at the start of the loop.

I chose to use the second option because I think the code is cleaner that way, though it is arguably less readable.

import numpy as np
from scipy.integrate import quad_vec


def f0(x, N):
    res = np.sin((N - 2)*np.pi*x)
    res[N == 1] = (1-x)
    res[N == 2] = x

    return res


def f1(x, N):
    res = np.cos((N - 2)*np.pi*x)*(N-2)*np.pi
    res[N == 1] = -1
    res[N == 2] = 1

    return res


def f00(x, M, N):
    return f0(x, M)*f0(x, N)


def f11(x, M, N):
    return f1(x, M)*f1(x, N)


def f10(x, M, N):
    return f1(x, M)*f0(x, N)


def f01(x, M, N):
    return f0(x, M)*f1(x, N)


NMHT_0 = 12

M, N = np.meshgrid(np.arange(0, NMHT_0), np.arange(0, NMHT_0))
I00, _ = quad_vec(f00, 0, 1, args=(M, N))
I11, _ = quad_vec(f11, 0, 1, args=(M, N))
I10, _ = quad_vec(f10, 0, 1, args=(M, N))
I01, _ = quad_vec(f01, 0, 1, args=(M, N))


a = 1
b = 1
A = np.random.rand(5, 5)

Kuu_s = np.zeros((NMHT_0**2, NMHT_0**2))

k = -1
for i in range(NMHT_0):
    for j in range(NMHT_0):
        k += 1
        q = -1
        for m in range(NMHT_0):
            for n in range(NMHT_0):
                q += 1
                Kuu_s[k, q] = (
                    A[1, 1]*I11[i, m]*I00[j, n]/a**2
                    + A[3, 3]*I00[i, m]*I11[j, n]/b**2
                    + A[1, 3]*(I10[i, m]*I01[j, n]
                               + I01[i, m]*I10[j, n])/a/b
                )

np.savetxt("matrixkuu2.txt", Kuu_s, delimiter="  ", fmt="%1.3e")

You'll also note that I got rid of the initial double for loops used for integration in favor of properly using quad_vec for evaluating the integrals. There might be a way to get rid of the large nested loop at the end, but I didn't bother. It will only be a problem (i.e. it will be slow) if NHMT_0 gets very large.