pcolormesh绘图维度不匹配问题求助（含有限差分代码）

问题描述 #

我用有限差分法求解方程组，输出的T_results和c_results数组末尾存在无效0值，于是对数组做切片处理得到：

T_results_cleaned = T_results[:, :-1]
c_results_cleaned = c_results[:, :-1]

但使用plt.pcolormesh绘制这些切片数组时，触发错误：

TypeError: Dimensions of C (50, 49) should be one smaller than X(50) and Y(50) while using shading='flat' see help(pcolormesh)

尝试将z和t数组也切片为z[:-1]、t[:-1]以匹配维度，问题仍未解决。

完整代码 #

import numpy as np
import matplotlib.pyplot as plt
import time 
import pandas as pd
nz = 50
dz = 1 / nz
nt = 50
t = np.linspace(0, 4, nt)

p = np.array([0.01, 0.5, 1.0])

# Initialize concentration and temperature arrays
C = np.zeros(nz)
T = np.zeros(nz)

# Set initial condition for concentration
C[0] = 1.0

# Initialize arrays to store results
c_results = np.zeros((nt, nz))
T_results = np.zeros((nt, nz))

# Define the Arrhenius rate constant function
def k(temperature):
    arrh_expr = 4.68 * np.exp(-806 / ((600 - 273) * temperature + 273))
    return arrh_expr

# Calculate time derivative at specific spatial points
def rhsC(c0, c1, c2, T1, dz, p):
    return (p[0] * (c2 - 2 * c1 + c0) / dz**2 - 
            p[1] * (c2 - c0) / (2 * dz) - 
            p[2] * k(T1) * c1**2)

def rhsT(T0, T1, T2, c1, dz, p):
    return (0.01 * (T2 - 2 * T1 + T0) / dz**2 - 
            0.5 * (T2 - T0) / (2 * dz) + 
            2.45 * k(T1) * c1**2)

# Time-stepping loop
for it in range(nt):
    c_results[it, :] = C
    T_results[it, :] = T

    T_results_cleaned = T_results[:, :-1] #select all rows and all columns except the last column
    c_results_cleaned = c_results[:, :-1] #select all rows and all columns except the last column
    
    for iz in range(1, nz - 1):
        # Calculate concentration and temperature derivatives
        dc_dt = rhsC(C[iz - 1], C[iz], C[iz + 1], T[iz], dz, p)
        dT_dt = rhsT(T[iz - 1], T[iz], T[iz + 1], C[iz], dz, p)
        
        # Update concentration and temperature using forward Euler
        C[iz] += dz * dc_dt
        T[iz] += dz * dT_dt

# Create a mesh for plotting
z = np.linspace(0, 1, nz)

# Plot concentration and temperature profiles
plt.pcolormesh(z, t, c_results_cleaned)
plt.colorbar(label='Concentration')
plt.xlabel('z')
plt.ylabel('t')
plt.title('Concentration Profile')
plt.show()

plt.pcolormesh(z, t, T_results_cleaned)
plt.colorbar(label='Temperature')
plt.xlabel('z')
plt.ylabel('t')
plt.title('Temperature Profile')
plt.show()

解决方案 #

1. 修正切片时机（先解决代码冗余问题） #

你当前在时间循环内部重复执行切片操作，完全没有必要，应该将切片代码移到时间循环结束之后，避免无效计算：

# 时间循环结束后再做切片
T_results_cleaned = T_results[:, :-1]
c_results_cleaned = c_results[:, :-1]

2. 解决pcolormesh维度匹配问题 #

plt.pcolormesh的shading='flat'模式（默认）要求颜色数组C的维度为(N-1, M-1)，其中X（z）的长度为M，Y（t）的长度为N（对应网格的顶点）。你的切片后数组是(50,49)，和X(50)、Y(50)不匹配，有两种可行方案：

方案A：切换到shading='nearest'模式

该模式允许C的维度与X、Y的点数量直接匹配（即C的行对应Y的点，列对应X的点），无需调整数组维度：

# 绘制浓度
plt.pcolormesh(z, t, c_results_cleaned, shading='nearest')
plt.colorbar(label='Concentration')
plt.xlabel('z')
plt.ylabel('t')
plt.title('Concentration Profile')
plt.show()

# 绘制温度
plt.pcolormesh(z, t, T_results_cleaned, shading='nearest')
plt.colorbar(label='Temperature')
plt.xlabel('z')
plt.ylabel('t')
plt.title('Temperature Profile')
plt.show()

方案B：调整网格与结果数组的维度匹配flat模式

如果坚持使用shading='flat'，需要让C的维度为(nt-1, nz-1)，同时将时间数组t保持原长度，结果数组去掉最后一个时间步的数据：

# 调整结果数组：去掉最后一个时间步的结果
c_results_flat = c_results_cleaned[:-1, :]
T_results_flat = T_results_cleaned[:-1, :]

# 绘制浓度
plt.pcolormesh(z, t, c_results_flat)
plt.colorbar(label='Concentration')
plt.xlabel('z')
plt.ylabel('t')
plt.title('Concentration Profile')
plt.show()

# 绘制温度
plt.pcolormesh(z, t, T_results_flat)
plt.colorbar(label='Temperature')
plt.xlabel('z')
plt.ylabel('t')
plt.title('Temperature Profile')
plt.show()

方案C：使用空间中点坐标匹配结果

如果你的切片后的结果对应空间网格的中点，可以生成中点坐标数组来匹配：

# 生成空间中点坐标
z_mid = (z[:-1] + z[1:]) / 2

# 绘制浓度
plt.pcolormesh(z_mid, t, c_results_cleaned, shading='nearest')
plt.colorbar(label='Concentration')
plt.xlabel('z')
plt.ylabel('t')
plt.title('Concentration Profile')
plt.show()

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