Matlab代码中矩阵p第11列及后续列全为0的问题排查

运行下方Matlab代码后，矩阵p的第11列及后续所有列变为0，且计算出的p值与预期的0.94左右或更低不符，无法定位问题原因。

clear all; clc; close all;

tic;

z1 = 1;
z1 = 2;
z2 = 2*z1;
z = [z1,z2];
la1 = 1/3;
la2 = 1/3;
la = [la1,la2];
z_ave = (z1*la2 + z2*la1)/(la1 + la2);
sigma = 0.75; 
theta = sigma/(sigma-1);
pe = 1.68;
xi = 0.01:0.01:0.2;

pe = linspace(1,1.5,100);
p = zeros(length(xi),length(pe))
for i = 1:20
    for j = 1:10
p(i,j) = (1+xi(i)*pe(j)^theta)^(1/theta)
end
end

plot(pe,p(1,:))

1. 列全为0的直接原因 #

代码内层循环j的范围被硬编码为1:10，但pe是长度为100的数组，p初始化时是20×100的零矩阵，这导致第11到100列从未被赋值，始终保持初始的0值。

修复方法：将内层循环范围改为1:length(pe)，同时外层循环用length(xi)替代硬编码的20，提升代码灵活性：

for i = 1:length(xi)
    for j = 1:length(pe)
        p(i,j) = (1+xi(i)*pe(j)^theta)^(1/theta);
    end
end

2. 数值与预期不符的原因 #

计算式中的theta为负数：sigma=0.75，因此theta=0.75/(0.75-1)=-3，1/theta=-1/3。代入后实际计算逻辑为：

p(i,j) = (1 + xi(i) / (pe(j)^3))^(-1/3) = 1 / (1 + xi(i)/(pe(j)^3))^(1/3)

当pe(j)在1到1.5区间时，结果始终接近1（例如pe=1.5、xi=0.2时，结果≈0.98），与预期的0.94左右不符。

建议：重新核对公式推导逻辑，比如是否theta的计算应为(sigma-1)/sigma，或是公式结构存在错误。若为CES需求函数，当替代弹性sigma<1时，需确认公式形式是否需要调整。

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