Matlab中矩阵与向量间的动能求和实现问题求助

我在大学项目中需要在Matlab中实现矩阵与向量间的动能求和，代码前半部分经教授确认正确，但在动能计算环节遇到问题。

动能公式：每个碎片的动能为 $E_{k,i} = \frac{1}{2} m_i | \boldsymbol{v} + \Delta \boldsymbol{V}{c,i} |^2$，总动能为所有碎片动能的总和，其中$\boldsymbol{v}$是初始速度向量，$\Delta \boldsymbol{V}{c,i}$是第i个碎片的修正速度向量，$| \cdot |$表示向量的模长。

原代码如下：

%% Popolazione di detriti

% Costanti
Mv=10000; % 10000 kg massa nominale velivolo
V_0 = 200*.3048; % velocità al distacco dalla WKT in m/s a 14325 m

N=1000
D=rand(N,1); % vettore debris di 100 numeri casuali tra 0 e 1
somma=sum(D); 
m_i=(D/somma)*Mv; %mass casuale dei detriti
ver=sum(m_i); % verifica che la somma dei dei 100 pezzi razionalizzati restituisce il peso del velivolo

vx_i = randn(N,1); % componenti random di velocità nelle tre direzioni
vy_i = randn(N,1); 
vz_i= randn(N,1);
DeltaV_iStar= [vx_i,vy_i,vz_i]; % matrix velocity i-esima debris
    %% Momentum
DeltaQ_err=zeros(1,3);
DeltaQ=zeros(N,3); %inizializzo matrice 
for k=1:N
    DeltaQ(k,:)=(m_i(k)*DeltaV_iStar(k,:));  
    DeltaQ_err=DeltaQ_err+DeltaQ(k,:);
end
DeltaQ_err


DeltaV_err = DeltaQ_err/Mv ; % errore da togliere agli incrementi iniziali
DeltaV_c = DeltaV_iStar-DeltaV_err;

DeltaQ_err2=zeros(1,3);
DeltaQ=zeros(N,3); %inizializzo matrice 
for k=1:N
  DeltaQ(k,:)=(m_i(k)*DeltaV_c(k,:));  
  DeltaQ_err2=DeltaQ_err2+DeltaQ(k,:);
end
DeltaQ_err2


%% Kinetic Energy
V_element= randn(3,1); % componenti della velocità iniziale
B = V_element/norm(V_element) 
v=V_0*B
vlength= norm(v);

prodotto=0;
Ek_d=0;
for k=1:N 
    prodotto(k,:)=.5*[m_i(k)*(v(k)+DeltaV_c(k,:)).^2];
    Ek_d=Ek_d+prodotto(k,:);
end
Ek_d

运行时报错：

Unable to perform assignment because the size of the left side is
1-by-1 and the size of the right side is 1-by-3.

Error in Debris_Footprint (line 59)
prodotto(k,:)=.5*[m_i(k)*(v(k)+DeltaV_c(k,:)).^2];

我尝试用for循环实现，但得到向量结果，而动能是标量，需要修正代码以正确计算标量形式的动能总和。

1. 索引错误：v是3×1的初始速度向量，循环中v(k)当k>3时会触发索引越界，因为v只有3个元素，正确的做法是用整个v向量和每个碎片的DeltaV_c(k,:)相加。
2. 向量运算错误：(v + DeltaV_c(k,:)).^2会对每个分量单独平方，得到1×3的向量，但动能需要的是速度向量模长的平方（即分量平方和），而非分量平方组成的向量。
3. 变量维度不匹配：prodotto(k,:)被初始化为标量0，后续赋值1×3的向量会导致维度不匹配报错。

方案1：修正for循环版本 #

%% Kinetic Energy
V_element= randn(3,1); % componenti della velocità iniziale
B = V_element/norm(V_element);
v=V_0*B;
vlength= norm(v);

Ek_d=0;
for k=1:N 
    % 计算第k个碎片的合速度向量
    v_total = v' + DeltaV_c(k,:); % v是3×1，转置为1×3和DeltaV_c(k,:)匹配
    % 计算合速度的模长平方，再乘以0.5*m_i(k)得到单个碎片动能，累加到总动能
    Ek_d = Ek_d + 0.5 * m_i(k) * dot(v_total, v_total);
end
Ek_d

方案2：向量化实现（更高效，无需循环） #

Matlab擅长向量化运算，直接用矩阵操作替代循环，速度更快：

%% Kinetic Energy
V_element= randn(3,1); % componenti della velocità iniziale
B = V_element/norm(V_element);
v=V_0*B;
vlength= norm(v);

% 将v转为1×3的行向量，与DeltaV_c（N×3）逐行相加，得到每个碎片的合速度矩阵
v_total_matrix = repmat(v', N, 1) + DeltaV_c;
% 计算每行的模长平方：sum(.*,2)表示按行求和
speed_squared = sum(v_total_matrix .* v_total_matrix, 2);
% 总动能：0.5 * 点积（质量向量和速度平方向量）
Ek_d = 0.5 * m_i' * speed_squared;
Ek_d

两种方案都能正确计算标量形式的总动能：

• 方案1通过循环逐个计算每个碎片的动能，逻辑直观，适合理解；
• 方案2利用Matlab的向量化特性，避免循环，计算效率更高，尤其当N很大时优势明显。

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