MATLAB 根函数在 MATLAB 和 Simulink 中的行为不同？

Question

我在 MATLAB 中实现了以下用户定义函数：

function Q = Calc_Q(Head, freq)
b6 = [3.7572E-07 -1.5707E-05 6.0490E-03 5.0018E-02 2.1180E-01];
b5 = [-9.0927E-06 8.9033E-04 -3.2415E-02 5.4525E-01 -8.1649E+00] / 10e2;
b4 = [7.5172E-06 -5.6565E-04 1.0024E-02 3.5888E-01 3.8894E-02] / 10e5;
b3 = [-4.8767E-06 4.8787E-04 -1.3311E-02 -1.2189E-01 -5.3522E+00] / 10e8;
b2 = [5.9227E-06 -8.1716E-04 3.5392E-02 -4.5413E-01 1.9547E+00] / 10e11;
b1 = [-2.0004E-06 2.9027E-04 -1.3754E-02 2.3490E-01 -1.2363E+00] / 10e14;

a = [polyval(b1,abs(freq)), polyval(b2, abs(freq)), polyval(b3, abs(freq)), polyval(b4, abs(freq)), polyval(b5, abs(freq)), polyval(b6, abs(freq)) - Head];

Q_roots = roots(a);
%Delete roots with imaginary part
i = 1;
while i <= length(Q_roots)
    if(imag(Q_roots(i)) ~= 0)
        Q_roots(i) = [];
        i = i - 1;
    end
    i = i + 1;
end
%Delete roots with real part greater then 3100
i = 1;
while i <= length(Q_roots)
    if(Q_roots(i) >= 3100 || Q_roots(i) < 0)
        Q_roots(i) = [];
        i = i - 1;
    end
    i = i +1;
end

if freq < 0
    Q = real(Q_roots(1)) * -1;
else
    Q = real(Q_roots(1));        
end
end

当我在 Matlab 中调用此函数时，它工作正常。但是，如果我在 simulink 中使用这个确切的代码作为 MATLAB 函数，它将停止工作。 （实际上它有效，但输出始终为零。）

我确实怀疑问题可能是什么。在调试模式下运行脚本时，我无法查看 Q_roots 的结果（它只是不显示任何内容）。

Q_roots = roots(a);

有什么想法吗？

Answer 1

问题很可能是由于您的逻辑消除了虚部中没有完全零的任何根。这是一种数学思维方式，在数值上并不能很好地发挥作用，至少在一般情况下是这样。可能在这两种情况下都找到了所有的根源（没有限制意味着其他情况），但在 Simulink 和代码生成中，问题被视为一个复杂的问题，并且一些根源可能会以微小的虚部返回。如果根的虚部不完全为零，则不要删除根，而是消除虚部在数值上无关紧要的根，或者相对于实部非常小相对或完全非常小.像

tol = 10*eps(class(Q_roots));
keepers = abs(imag(Q_roots)) < tol*max(abs(real(Q_roots)),1) & ...
    real(Q_roots) >= 0 & real(Q_roots) <= 3100;
Q_roots = Q_roots(keepers);

会一举解决所有删除问题。我在这里使用 10*eps 作为容差。

但如果您只需要第一个符合条件的根，那么您可以这样做：

Q = nan('like',a);
tol = 10*eps(class(a));
for k = 1:numel(Q_roots)
    r = real(Q_roots(k));
    if abs(imag(Q_roots(k))) < tol*max(abs(r),1) && r >= 0 && r <= 3100;
        Q = r;
        break
    end
end

if freq < 0
    Q = -Q;
end

Answer 2

好的，我发现了问题。

来自不同的论坛：

嗨，科斯敏，

我查看了嵌入式 MATLAB 根的实现 功能块 (\toolbox\eml\lib\matlab\polyfun\roots.m)。 那里有说明：

% 限制：% 输出总是可变大小。 % 输出是 总是很复杂。 % Roots 的顺序可能与 MATLAB 不同。 % 条件差的多项式的根可能与 MATLAB 不匹配。最后 这句话让你头疼（是的，你的多项式是 条件很差！）。如果你看一下你会看到的情节，那 曲线几乎不接触 x 轴。

不过我有一个建议：值 -z/b 是（非常）好 您正在寻找的根的近似值...？

提图斯

http://www.mathworks.com/matlabcentral/answers/25624-roots-in-simulink

显然，simulink 中的根函数并不总是找到给定多项式的所有根。 这是不幸的，也不容易解决。不过我确实找到了解决方案。

对于我必须解决的所有不同多项式，我知道我感兴趣的根的区间（[-3000, 3000]）。

我基本上只是从-3000到3000走50步，直到函数降到0以下。然后我知道了根的近似解。我将此近似值用作 Newton-Raphson 方法的种子。

对于我必须求解的所有多项式，使用给定种子直接实现牛顿拉夫森方法不起作用，因为有时它会迭代到不同的根（我不感兴趣的根）。

代码如下：

function Q = Calc_Q(Head, freq)
    b6 = [3.7572E-07 -1.5707E-05 6.0490E-03 5.0018E-02 2.1180E-01];
    b5 = [-9.0927E-06 8.9033E-04 -3.2415E-02 5.4525E-01 -8.1649E+00] / 10e2;
    b4 = [7.5172E-06 -5.6565E-04 1.0024E-02 3.5888E-01 3.8894E-02] / 10e5;
    b3 = [-4.8767E-06 4.8787E-04 -1.3311E-02 -1.2189E-01 -5.3522E+00] / 10e8;
    b2 = [5.9227E-06 -8.1716E-04 3.5392E-02 -4.5413E-01 1.9547E+00] / 10e11;
    b1 = [-2.0004E-06 2.9027E-04 -1.3754E-02 2.3490E-01 -1.2363E+00] / 10e14;

    %coeff for the polynominal
    a = [polyval(b1,abs(freq)), polyval(b2, abs(freq)), polyval(b3, abs(freq)), polyval(b4, abs(freq)), polyval(b5, abs(freq)), polyval(b6, abs(freq)) - Head];

    %coeff for the derrivative of polynominal
    da = [5*a(1) 4*a(2) 3*a(3) 2*a(4) a(5)];

    Q = -3000;
    %Search for point where function goes below 0
    while (polyval(a, Q) > 0)
        Q = Q + 25;
    end    
    error_max = 0.01
    iter_counter = 1;
    while abs(polyval(a,Q)) >= error_max && iter_counter <= 1000
        Q = Q - polyval(a, Q)/polyval(da, Q);
        iter_counter = iter_counter + 1;
        error = abs(polyval(a,Q));
    end
    if(freq < 0)
        Q = Q * - 1;
    end
end