在Matlab中获得使积分等于零的常数答案

【问题标题】：Obtaining the constant that makes the integral equal to zero in Matlab在Matlab中获得使积分等于零的常数
【发布时间】：2016-03-17 01:18:28
【问题描述】：

我正在尝试编写一个 MATLAB 程序，我已经到了需要执行以下操作的地步。我有这个等式：

我必须找到常数“Xcp”的值（大于零），也就是使积分等于零的值。

为了做到这一点，我编写了一个循环，其中 Xcp 的值在每次迭代中以小增量前进，并执行积分并检查它是否为零，如果它达到零，则循环结束并且 Xcp 为与此值一起存储。

但是，我认为这不是完成这项任务的有效方法。运行时间增加了很多，因为这个循环很长，每次都要进行积分，积分限制替换。

在 Matlab 中是否有更聪明的方法来获得更好的代码效率？

P.S.：我使用conv() 将两个多项式相乘。因为 cl(x) 和 (x-Xcp) 都是多项式。

编辑：一段代码。

p = [1 -Xcp];    % polynomial (x-Xcp)
Xcp=0.001;
i=1;
found=false;
while(i<=x_te && found~=true)      % Xcp is upper bounded by x_te
   int_cl_p = polyint(conv(cl,p));
   Cm_cp=(-1/c^2)*diff(polyval(int_cl_p,[x_le,x_te]));
   if(Cm_cp==0)
     found=true;
   else
     Xcp=Xcp+0.001;
   end
end

这是我用来运行本节的代码。另一个问题是我必须针对不同的情况（不同的 cl 函数）这样做，因此代码更慢。

【问题讨论】：

最好的方法是借助矢量化来做到这一点。出示您的代码，我们会尽力帮助您！
您可以重新排列并从积分中取出 Xcp 以获得 Xcp = A/B 其中 A = Integral( x*cl(x)) 和 B = Integral(cl(x))
这是个好主意，谢谢@dmuir

标签： matlab math integration

【解决方案1】：

据我所知，您需要求解 X_CP 的方程。我建议为此使用符号求解器。对于大型多项式，这不是最有效的方法，但对于 20 次多项式，它需要不到 1 秒的时间。我并不是说这个解决方案是最快的，但它提供了通用解决方案。如果你的多项式每次迭代都不会改变，那么你可以多次使用这个通用解决方案，而不用花时间计算积分。

因此，就xLE 和xTE 而言，通用符号解决方案是使用以下方法获得的：

syms xLE xTE c x xCP
a = 1:20;
%//arbitrary polynomial of degree 20
cl = sum(x.^a.*randi([-100,100],1,20));
tic
eqn = -1/c^2 * int(cl * (x-xCP), x, xLE, xTE) == 0;
xCP = solve(eqn,xCP);
pretty(xCP)
toc

Elapsed time is 0.550371 seconds.

您可以进一步使用matlabFunction 来寻找数值解：

xCP_numerical = matlabFunction(xCP);
%// we then just plug xLE = 10 and xTE = 20 values into function
answer = xCP_numerical(10,20)
answer =  
        19.8038

对代码的轻微修改可以让您将其用于通用系数。

希望有帮助

【讨论】：

+1，一个很好的符号工具箱应用。我建议使用 subs 来评估它，而不是 matlabFunction。这样您就不会失去精度。
非常感谢。我会这样做，因为它变得更有效率，而且我的多项式是 9 次的，所以它会更快。感谢您的帮助。
对多项式积分使用符号求解器是一种耻辱。当方程是线性的时，使用另一个求解器来获得未知数是另一种耻辱。
@YvesDaoust 如果您很聪明，请提出您的解决方案。将您居高临下的 cmets 留在此线程之外
关注@Steve 并使用 int[a,b] x^kdx=(b^(k+1)-a^(k+1))/(k+1)。你得到一个封闭的公式，不需要软件。

【解决方案2】：

如果乘以-1/c^2，则可以重新排列为

并根据您的喜好进行集成。由于c_l 是多项式阶N，如果它是在MATLAB 中使用polyval 的常用符号定义的，其中系数存储在向量a 中，这样

那么集成很简单：

MATLAB 代码可能看起来像这样

int_cl_p = polyint(cl);
int_cl_x_p = polyint([cl 0]);
X_CP = diff(polyval(int_cl_x_p,[x_le,x_te]))/diff(polyval(int_cl_p,[x_le,x_te]));

【讨论】：