You need to enable JavaScript to run this app.
优惠活动
大模型
产品
解决方案
定价
更多

MATLAB浓度数组应为840×1却返回840×840的问题求助

avatar
API 探路者
2026-7-29
多湖泊污染传播模拟MATLAB代码数组维度异常排查

我正在运行一段模拟多湖泊不同污染质量传播的MATLAB代码，计算各湖泊污染浓度后尝试绘图，但发现浓度数组应为840×1，实际却返回840×840。以下是我的代码：

% SECTION 1 - Load variables
clear
close all
clc
% Replace the path within quotation marks below to the folder where the
% data is stored in your computer
cd '/Users/brook/Project_Support_Files/Data/Flow Data/'
FL = dir('*.csv');
% Read the spreadsheet with dates and transfom it into an array
t = readtable('Date.xlsx');
t = table2array(t);
% surface areas and volumes

A1 = 8.1925e10;                     % Superior
V01 = 12000*1E9;

A2 = 1.1685e11;                     % Mi-Huron
V02 = (3500+4900)*1E9;

A3 = 2.5404e10;                     % Erie
V03 = 480 * 1E9;

A4 = 1.9121e10;                     % Ontario
V04 = 1640 * 1E9;
% load diversions
% Terms multiplied convert data from m3/s to mm over the lakes surface area
D1 = readmatrix('D1.csv')*1000*60*60*24*30/A1;
D2 = readmatrix('D2.csv')*1000*60*60*24*30/A2;
D3 = readmatrix('D3.csv')*1000*60*60*24*30/A3;
% load connecting flows
F1 = readmatrix('Fo1.csv')*1000*60*60*24*30/A1;
F2 = readmatrix('Fo2.csv')*1000*60*60*24*30/A2;
F3 = readmatrix('Fo3.csv')*1000*60*60*24*30/A3;
F4 = readmatrix('Fo4.csv')*1000*60*60*24*30/A4;
% load evaporation
E1 = readmatrix('E1.csv');
E2 = readmatrix('E2.csv');
E3 = readmatrix('E3.csv');
E4 = readmatrix('E4.csv');
% load precipitation
P1 = readmatrix('P1.csv');
P2 = readmatrix('P2.csv');
P3 = readmatrix('P3.csv');
P4 = readmatrix('P4.csv');
% load runoff
R1 = readmatrix('R1.csv');
R2 = readmatrix('R2.csv');
R3 = readmatrix('R3.csv');
R4 = readmatrix('R4.csv');
% compute ODEs 5.1 to 5.4
% Multiplication by A_i and 1e-3 converts to m3 per month
dV1 = (  R1 + P1 + D1 - F1 - E1) * A1                       * 1E-3;
dV2 = ( (R2 + P2 - D2 - F2 - E2) * A2 +       (F1 * A1) )   * 1E-3;
dV3 = ( (R3 + P3 - D3 - F3 - E3) * A3 +       (F2 * A2) )   * 1E-3;
dV4 = ( (R4 + P4      - F4 - E4) * A4 + (D3 + F3) * A3  )   * 1E-3;
% Integration loops. Timestep = 1month
% Notice that V_i (in this case) is simply the cumulative sum of the
% differentials. For concentration solutions a Euler scheme needs to be
% implemented
V1 = zeros(size(dV1));
V2 = zeros(size(dV2));
V3 = zeros(size(dV3));
V4 = zeros(size(dV4));
for in = 1 : length(dV1)
    V1(in) = sum(dV1(1:in));
    V2(in) = sum(dV2(1:in));
    V3(in) = sum(dV3(1:in));
    V4(in) = sum(dV4(1:in));
end
clear in
V1 = V1 + V01;
V2 = V2 + V02;
V3 = V3 + V03;
V4 = V4 + V04;
%% Plot volume results
% Use this template of plots for your results for concentrations
% Notice that the subplot function helps you plot several graphs in a
% single figure
figure()
subplot(2,2,1)
plot(t,V1*1e-9,'k', 'linewidth', 2)
ylabel('V [km^3]')
title('Superior')
grid on
subplot(2,2,2)
plot(t,V2*1e-9,'b', 'linewidth', 2)
title('Michigan-Huron')
grid on
subplot(2,2,3)
plot(t,V3*1e-9,'c', 'linewidth', 2)
xlabel('Year')
title('Erie')
ylabel(' V [km^3]')
grid on
subplot(2,2,4)
plot(t,V4*1e-9,'r', 'linewidth', 2)
xlabel('Year')
title('Ontario')
grid on
% Concentration of run-off pollution
CR1=0;
CR2=0;
CR3=0;
CR4=0;
%Converting R to litre per month from mm over surface area
R1=R1*A1*(10^(-3))*1000;
R2=R2*A2*(10^(-3))*1000;
R3=R3*A3*(10^(-3))*1000;
R4=R4*A4*(10^(-3))*1000;
%Converting F to litre per month from mm over surface area
F1=F1*A1*(10^(-3))*1000;
F2=F2*A2*(10^(-3))*1000;
F3=F3*A3*(10^(-3))*1000;
F4=F4*A4*(10^(-3))*1000;
%Converting V to litres from m3
V1=V1*1000;
V2=V2*1000;
V3=V3*1000;
V4=V4*1000;
% M values in micrograms (for initial mass of pollution 75000 tonnes )
M01=75000*(10^6)*(10^6);
M02=0;
M03=0;
M04=0;
%Time step
h=1;
t=1:h:length(t);
M1=zeros(size(t));
M2=zeros(size(t));
M3=zeros(size(t));
M4=zeros(size(t));
M1(1)=M01;
M2(1)=M02;
M3(1)=M03;
M4(1)=M04;
% Approximate solution for next values of M1,M2, M3 and M4:
for i=1:(length(t)-1)
dM1=CR1*R1-(M1(i)./V1).*F1;
M1(i+1)=M1(i)+dM1(i)*h;
dM2=CR2*R2+(M1(i)./V1).*F1-(M2(i)./V2).*F2;
M2(i+1)=M2(i)+dM2(i)*h;
dM3=CR3*R3+(M2(i)./V2).*F2-(M3(i)./V3).*F3;
M3(i+1)=M3(i)+dM3(i)*h;
dM4=CR4*R4+(M3(i)./V3).*F3-(M4(i)./V4).*F4;
M4(i+1)=M4(i)+dM4(i)*h;
end
C1=(M1./V1);

C2=(M2./V2);

C3=(M3./V3);

C4=(M4./V4);
% M values in micrograms ( for initial mass of pollutant 50,000 tonnes)
M012=50000*(10^6)*(10^6);
M022=0;
M032=0;
M042=0;
M12=zeros(size(t));
M22=zeros(size(t));
M32=zeros(size(t));
M42=zeros(size(t));
M12(1)=M012;
M22(1)=M022;
M32(1)=M032;
M42(1)=M042;
% Approximate solution for next values of M12, M22, M32 and M42
for i=1:(length(t)-1)
dM12=CR1*R1-(M12(i)./V1).*F1;
M12(i+1)=M12(i)+dM12(i)*h;
dM22=CR2*R2+(M12(i)./V1).*F1-(M22(i)./V2).*F2;
M22(i+1)=M22(i)+dM22(i)*h;
dM32=CR3*R3+(M22(i)./V2).*F2-(M32(i)./V3).*F3;
M32(i+1)=M32(i)+dM32(i)*h;
dM42=CR4*R4+(M32(i)./V3).*F3-(M42(i)./V4).*F4;
M42(i+1)=M42(i)+dM42(i)*h;
end
C12=(M12./V1);
C22=(M22./V2);
C32=(M32./V3);
C42=(M42./V4);
% M values in micrograms (for initial mass pollution 25000 tonnes)
M013=25000*(10^6)*(10^6);
M023=0;
M033=0;
M043=0;
M13=zeros(size(t));
M23=zeros(size(t));
M33=zeros(size(t));
M43=zeros(size(t));
M13(1)=M013;
M23(1)=M023;
M33(1)=M033;
M43(1)=M043;
% Approximate solution for next values of M13, M23, M33 and M43
for i=1:(length(t)-1)
dM13=CR1*R1-(M13(i)./V1).*F1;
M13(i+1)=M13(i)+dM13(i)*h;
dM23=CR2*R2+(M13(i)./V1).*F1-(M23(i)./V2).*F2;
dM23=reshape(dM23',[840,1]);
M23(i+1)=M23(i)+dM23(i)*h;
dM33=CR3*R3+(M23(i)./V2).*F2-(M33(i)./V3).*F3;
dM33=reshape(dM33',[840,1]);
M33(i+1)=M33(i)+dM33(i)*h;
dM43=CR4*R4+(M33(i)./V3).*F3-(M43(i)./V4).*F4;
dM43=reshape(dM43',[840,1]);
M43(i+1)=M43(i)+dM43(i)*h;
end
C13=(M13./V1);
C23=(M23./V2);
C33=(M33./V3);
C43=(M43./V4);
t = readtable('Date.xlsx');
t = table2array(t);
figure()
subplot(2,2,1)
plot(t,C1,'k',t,C12,'y',t,C13, 'c','linewidth',2)
ylabel('Concentration [micrograms per litre]')
title('Superior')
grid on
subplot(2,2,2)
plot(t,C2,'k',t,C22,'y',t,C23, 'c','linewidth',2)
title('Michigan-Huron')
grid on
subplot(2,2,3)
plot(t,C3,'k', t,C32,'y', t,C33, 'c','linewidth',2)
xlabel('Year')
title('Erie')
ylabel(' concentration [micrograms per litre]')
grid on
subplot(2,2,4)
plot(t,C4,'k', t,C42,'y',t,C43,'c','linewidth',2)
xlabel('Year')
title('Ontario')
grid on

维度异常原因排查

  • 核心问题：广播机制导致的矩阵生成
    代码中M1(i)是标量，V1为840×1向量，(M1(i)./V1)得到840×1向量；当它与同样是840×1的F1执行.*运算时，MATLAB的广播机制会将两个列向量分别转置后进行外积运算，最终生成840×840的矩阵，而非预期的逐元素乘积（840×1向量）。后续dM2、dM3等计算重复此错误，导致M数组被逐步填充为矩阵，最终C = M./V也变成840×840矩阵。

  • 具体出错示例
    以第一个循环中的dM1计算为例：

    dM1=CR1*R1-(M1(i)./V1).*F1;

    (M1(i)./V1)是840×1，F1是840×1，两者的.*运算触发广播生成矩阵，后续取dM1(i)只是提取矩阵的第i个元素，但dM1本身已经是矩阵，导致M1数组逐渐变为矩阵。

  • 解决方案

    1. 统一向量维度：检查加载的F1、R1等数据，若为行向量，转置为列向量：
      F1 = readmatrix('Fo1.csv')*1000*60*60*24*30/A1;
F1 = F1.'; % 转置为列向量
    2. 避免不必要的广播：确保所有元素-wise运算的操作数维度完全匹配，比如确认V1、F1等均为列向量，这样(M1(i)./V1).*F1会直接生成840×1的向量，而非矩阵。
    3. 移除无效的reshape操作：第三个循环中的reshape只是临时修复，解决根源问题后无需保留。

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

相关产品推荐
方舟 Agent Plan

超全模态模型 × Harness 升级，最新支持 Deepseek-V4.1-Flash、GLM-5.3 系列、Doubao-Seedream-5.0-pro、Kimi-K3 (部分)， 限时 9.9 元起

产品介绍相关文档
相关产品推荐
方舟 Agent Plan

超全模态模型 × Harness 升级，最新支持 Deepseek-V4.1-Flash、GLM-5.3 系列、Doubao-Seedream-5.0-pro、Kimi-K3 (部分)， 限时 9.9 元起

产品介绍相关文档
火山引擎 最新活动
全面开放
Doubao-Seedance-2.5全新模型升级

模型再升级,30秒超长叙事, 模态参考扩容

HOT
方舟Codingplan

模型自由，工具不限，最新支持 Deepseek-V4 系列、GLM-5.3 系列

HOT
方舟Agent plan

超全模态模型 × Harness 升级，最新支持 Deepseek-V4.1-Flash、GLM-5.3 系列

HOT
新用户特惠专场

大模型19元起，Al应用9.9元畅享，新人首购爆款尽享优惠

最近更新时间：2026.07.29 05:57:03