Python中基于[i][j]张量元素积分生成矩阵结果错误的求助

核心错误点 #

1. 被积函数参数错误：你代码里的被积函数用了(n-2)πx和(m-2)πx，但Mathcad中的表达式是(m-1)πx与(n-1)πx的正弦乘积，这是导致结果不正确的根本原因。
2. 循环范围与索引不匹配：
   • 你设置的循环是range(1, NMHT_0)，也就是m、n仅从1遍历到11（Python的range是左闭右开），漏掉了m=12和n=12的情况，导致矩阵最后一行/列始终为0。
   • Python数组是0-based索引，但你直接用m、n作为索引赋值，会导致第0行、第0列完全未被赋值，始终为0，和Mathcad的1-based索引逻辑不对应。

修正后的代码 #

方法1：修正数值积分逻辑 #

import numpy as np
import scipy.integrate as spi

NMHT_0 = 12
I00 = np.zeros((NMHT_0, NMHT_0))

# 修正被积函数的参数为(m-1)和(n-1)
def integrand(x, m, n):
    return np.sin((n - 1) * np.pi * x) * np.sin((m - 1) * np.pi * x)

# 遍历1到12的m、n，转换为0-based索引赋值
for m in range(1, NMHT_0 + 1):
    for n in range(1, NMHT_0 + 1):
        result, error = spi.quad(integrand, 0, 1, args=(m, n))
        I00[m-1, n-1] = result

# 打印结果
for row in I00:
    print(' '.join(f'{val:.3f}' for val in row))

方法2：利用正交性直接生成矩阵（更高效） #

由于三角函数的正交性，这个积分有解析解：

• 当m == n时，积分结果为0.5
• 当m != n时，积分结果为0

可以直接生成矩阵，无需数值积分：

import numpy as np

NMHT_0 = 12
# 创建对角元为0.5的矩阵
I00 = np.eye(NMHT_0) * 0.5

# 打印结果
for row in I00:
    print(' '.join(f'{val:.3f}' for val in row))

验证说明 #

两种方法生成的矩阵完全一致，方法2避免了数值积分的计算开销，适合大维度矩阵的生成。

内容的提问来源于stack exchange，提问作者Ayoub