Python用Matplotlib绘循环变量图遇错，需观测数值方法收敛性

报错原因 #

从报错信息来看，问题出在fi_1 = A[0]*sin(...)这一行：A是numpy矩阵（np.matrix）类型，A[0]返回的是一个(1,4)的二维矩阵，而sin(...)返回的是标量，二者相乘时维度不兼容，触发operands could not be broadcast together错误。

修正方案 #

有两种简单的解决方式：

方式1：将矩阵A转为一维数组 #

在计算完A后，将其转换为普通的numpy一维数组，这样索引就能直接得到标量：

A = np.dot(inverse_matrix_A, K_matrix)
# 转换为一维数组
A = np.asarray(A).ravel()

方式2：正确索引矩阵中的标量元素 #

因为np.matrix是二维结构，每个元素需要用二维索引提取标量，比如A[0, 0]、A[1, 0]等：

fi_1 = A[0,0]*sin(1*(W*i+B))
fi_2 = A[1,0]*sin(2*(W*i+B))
fi_3 = A[2,0]*sin(3*(W*i+B))
fi_4 = A[3,0]*sin(4*(W*i+B))

完整修正后代码（含绘图） #

以下是整合了修正方案（方式1）并添加Matplotlib绘图的完整代码：

1. 导入依赖库 #

import numpy as np
import matplotlib.pyplot as plt
from math import exp, pi, cos, sin, sqrt, log

2. 矩阵参数计算 #

x_a = 1  # 补充缺失的x_a定义，可按需修改
x_b = exp(pi/sqrt(3))

def C(x,y):
    C1 = cos(sqrt(3)/2 * log(x))/sqrt(x)
    C2 = cos(sqrt(3)/2 * log(y))/sqrt(y)
    C3 = sin(sqrt(3)/2 * log(x))/sqrt(x)
    C4 = sin(sqrt(3)/2 * log(y))/sqrt(y)
    return np.matrix('{} {};{} {}'.format(C1, C3, C2, C4))

matrix = C(x_a, x_b)
Y_a = 1
Y_b = 2

def Y(x,y):
    Y1 = x
    Y2 = y
    return np.matrix('{};{}'.format(Y1, Y2))

Y = Y(Y_a, Y_b)
inverse_array = np.linalg.inv(matrix)
C_mat = np.dot(inverse_array, Y)  # 重命名避免覆盖函数C
C1 = 1
C2 = 2 * exp(pi/(2*sqrt(3)))
k = 1/(exp(pi/sqrt(3))-1)
print("k:")
print(k)
b = 1 - k
W = pi/(exp(pi/sqrt(3))-1)
B = pi - W

3. 矩阵变量计算 #

X_1 = x_a + (x_b - x_a)/5
X_2 = X_1 + (x_b - x_a)/5
X_3 = X_2 + (x_b - x_a)/5
X_4 = X_3 + (x_b - x_a)/5
X_5 = X_4 + (x_b - x_a)/5

A11 = sin(1 * (W * X_1 + B)) * (1 / (X_1 * X_1) - 1 * 1 * W * W)
A12 = sin(1 * (W * X_2 + B)) * (1 / (X_2 * X_2) - 1 * 1 * W * W)
A13 = sin(1 * (W * X_3 + B)) * (1 / (X_3 * X_3) - 1 * 1 * W * W)
A14 = sin(1 * (W * X_4 + B)) * (1 / (X_4 * X_4) - 1 * 1 * W * W)

A21 = sin(2 * (W * X_1 + B)) * (1 / (X_1 * X_1) - 2 * 2 * W * W)
A22 = sin(2 * (W * X_2 + B)) * (1 / (X_2 * X_2) - 2 * 2 * W * W)
A23 = sin(2 * (W * X_3 + B)) * (1 / (X_3 * X_3) - 2 * 2 * W * W)
A24 = sin(2 * (W * X_4 + B)) * (1 / (X_4 * X_4) - 2 * 2 * W * W)

A31 = sin(3 * (W * X_1 + B)) * (1 / (X_1 * X_1) - 3 * 3 * W * W)
A32 = sin(3 * (W * X_2 + B)) * (1 / (X_2 * X_2) - 3 * 3 * W * W)
A33 = sin(3 * (W * X_3 + B)) * (1 / (X_3 * X_3) - 3 * 3 * W * W)
A34 = sin(3 * (W * X_4 + B)) * (1 / (X_4 * X_4) - 3 * 3 * W * W)

A41 = sin(4 * (W * X_1 + B)) * (1 / (X_1 * X_1) - 4 * 4 * W * W)
A42 = sin(4 * (W * X_2 + B)) * (1 / (X_2 * X_2) - 4 * 4 * W * W)
A43 = sin(4 * (W * X_3 + B)) * (1 / (X_3 * X_3) - 4 * 4 * W * W)
A44 = sin(4 * (W * X_4 + B)) * (1 / (X_4 * X_4) - 4 * 4 * W * W)

A_matrix = np.matrix('{} {} {} {};{} {} {} {};{} {} {} {};{} {} {} {}'.format(A11, A21, A31, A41, A12, A22, A32, A44, A13, A23, A33, A43, A14, A24, A34, A44))
inverse_matrix_A = np.linalg.inv(A_matrix)

def K(x, y, z):
    k_val = -(x/y + z/(y*y))
    return k_val

K_matrix = np.matrix('{};{};{};{}'.format(K(k, X_1, b), K(k, X_2, b), K(k, X_3, b), K(k, X_4, b)))
A = np.dot(inverse_matrix_A, K_matrix)
# 转换为一维数组解决维度问题
A = np.asarray(A).ravel()

vals = np.arange(x_a, x_b, 0.1)

4. 变量计算与绘图 #

kol, tochnoe, Eps = [], [], []
for i in vals:
    fi_0 = k*i+b
    fi_1 = A[0]*sin(1*(W*i+B))
    fi_2 = A[1]*sin(2*(W*i+B))
    fi_3 = A[2]*sin(3*(W*i+B))
    fi_4 = A[3]*sin(4*(W*i+B))
    kol.append(fi_0 + fi_1 + fi_2 + fi_3 + fi_4)
    tochnoe.append(C1 * cos(sqrt(3) / 2 * log(i))/sqrt(i) + C2 * sin(sqrt(3) / 2 * log(i))/sqrt(i))
    Eps.append(tochnoe[-1] - kol[-1])

# 绘制收敛性分析图
plt.figure(figsize=(12, 8))

# 数值解与精确解对比
plt.subplot(2,1,1)
plt.plot(vals, kol, label='数值解')
plt.plot(vals, tochnoe, label='精确解', linestyle='--')
plt.xlabel('x')
plt.ylabel('y')
plt.title('数值解与精确解对比')
plt.legend()

# 误差曲线（观测收敛性）
plt.subplot(2,1,2)
plt.plot(vals, Eps, label='误差', color='red')
plt.xlabel('x')
plt.ylabel('误差')
plt.title('数值方法收敛性误差')
plt.legend()

plt.tight_layout()
plt.show()

额外说明 #

• 补充了用户缺失的x_a定义（默认设为1，可根据实际需求修改）
• 将变量C重命名为C_mat，避免覆盖之前定义的函数C
• 添加了完整的Matplotlib绘图代码，直接运行即可观测收敛性

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