Matlab中变上限积分计算升力剪力与弯矩是否正确？

你的积分计算存在关键问题，导致结果不符合预期：

1. trapz函数的维度参数错误 #

你在调用trapz时传入了第三个参数2，但x(i:end)和pz(i:end)都是一维行向量，维度2指向列方向，而这两个向量的列长度仅为1，会直接导致积分计算出错。正确做法是去掉第三个参数（一维数组默认按行计算积分），或者显式传入1作为维度参数。

2. 物理趋势验证（辅助确认） #

从航空力学逻辑来看，剪力和弯矩的变化应该符合：

• 当x=0（翼根），剪力为整个半展升力的积分值，弯矩达到最大值；
• 当x=b/2（翼尖），剪力和弯矩均为0。

修正后的代码 #

% 需提前定义参数b和p0，示例值：b=2, p0=1000
b = 2;
p0 = 1000;

% 生成展向距离向量
x = linspace(0, b/2, 100);

% 计算升力分布pz(x)
pz = p0 * sqrt(1 - (x / (b/2)).^2);

% 预分配剪力和弯矩数组
Vz_lift = zeros(size(x));
My_lift = zeros(size(x));

% 循环计算每个位置的剪力和弯矩
for i = 1:length(x)
    % 修正trapz的维度参数（移除错误的2）
    Vz_lift(i) = trapz(x(i:end), pz(i:end));
    % 弯矩被积函数逻辑正确，同步修正trapz参数
    My_lift(i) = -trapz(x(i:end), (x(i:end) - x(i)) .* pz(i:end));
end

% 绘制剪力分布图
subplot(3,1,2)
plot(x, Vz_lift);
xlabel('展向距离 (x)');
ylabel('剪力分布 (N)');
title('半展升力产生的剪力分布');

% 绘制弯矩分布图
subplot(3,1,3)
plot(x, My_lift);
xlabel('展向距离 (x)');
ylabel('弯矩分布 (N·m)');
title('半展升力产生的弯矩分布');

解析解验证（用于对比数值结果） #

你可以用解析解验证修正后的数值结果是否正确：

• 剪力解析表达式：
  $$ V_z(x) = \frac{p_0 b}{4} \left( \frac{\pi}{2} - \arcsin\left(\frac{2x}{b}\right) - \frac{2x}{b}\sqrt{1-\left(\frac{2x}{b}\right)^2} \right) $$
• 弯矩解析表达式：
  $$ M_y(x) = -\frac{p_0 b^2}{24} \left( 3\left( \arcsin\left(\frac{2x}{b}\right) + \frac{2x}{b}\sqrt{1-\left(\frac{2x}{b}\right)^2} - \frac{\pi}{2} \right) + \left(1-\left(\frac{2x}{b}\right)^2\right){\frac{3}{2}} \right) $$

将解析解计算结果与修正后的数值积分结果对比，若趋势和数值一致，说明代码已修复。

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