使用 MuPAD 查找非线性方程组的所有解

Question

我的问题是是否有在 Matlab 脚本中使用 MuPAD 函数的好方法。背景是我有一个问题，我需要找到一组非线性方程的所有解。之前的解决方案是在 Matlab 中使用solve，它适用于我的一些模拟（即，一些输入集T），但并非总是如此。因此，我通过以下方式使用 MuPAD：

function ut1 = testMupadSolver(T)
% # Input T should be a vector of 15 elements

mupadCommand = ['numeric::polysysroots({' eq1(T) ' = 0,' ...
    eq2(T) '= 0},[u, v])'];
allSolutions = evalin(symengine, mupadCommand);
ut1 = allSolutions;

end

function strEq = eq1(T)
sT = @(x) ['(' num2str(T(x)) ')'];

strEq = [ '-' sT(13) '*u^4 + (4*' sT(15) '-2*' sT(10) '-' sT(11) '*v)*u^3 + (3*' ...
    sT(13) '-3*' sT(6) '+v*(3*' sT(14) '-2*' sT(7) ')-' sT(8) '*v^2)*u^2 + (2*' ...
    sT(10) '-4*' sT(1) '+v*(2*' sT(11) '-3*' sT(2) ')+v^2*(2*' sT(12) ' - 2*' ...
    sT(3) ')-' sT(4) '*v^3)*u + v*(' sT(7) '+' sT(8) '*v+' sT(9) '*v^2)+' sT(6)];

end

function strEq = eq2(T)
sT = @(x) ['(' num2str(T(x)) ')'];
strEq = ['(' sT(14) '-' sT(13) '*v)*u^3 + u^2*' '(' sT(11) '+(2*' sT(12) '-2*' sT(10) ...
    ')*v-' sT(11) '*v^2) + u*(' sT(7) '+v*(2*' sT(8) '-3*' sT(6) ')+v^2*(3*' sT(9) ...
    '-2*' sT(7) ') - ' sT(8) '*v^3) + v*(2*' sT(3) '-4*' sT(1) '+v*(3*' sT(4) ...
    '-3*' sT(2) ')+v^2*(4*' sT(5) ' - 2*' sT(3) ')-' sT(4) '*v^3)+' sT(2)];

end

我有两个疑问：

1) 为了使用 MuPAD，我需要将方程系统的两个方程重写为字符串，如上所示。有没有更好的方法来做到这一点，最好没有字符串步骤？

2) 关于格式输出；当

T = [0 0 0 0 0 0 0 0 0 0 1 0 1 0 1];

输出是：

testMupadSolver(T)

ans =

matrix([[u], [v]]) in {matrix([[4.4780323328249527319374854327354], [0.21316518769990291263811232040432]]), matrix([[- 0.31088044854742790561428736573347 - 0.67937835289645431373983117422178*i], [1.1103383836576028262792542770062 + 0.39498445715599777249947213893789*i]]), matrix([[- 0.31088044854742790561428736573347 + 0.67937835289645431373983117422178*i], [1.1103383836576028262792542770062 - 0.39498445715599777249947213893789*i]]), matrix([[0.47897094942962218512261248590261], [-1.26776233072168360314707025141]]), matrix([[-0.83524238515971910583152318717102], [-0.66607962429342496204955062300669]])} union solvelib::VectorImageSet(matrix([[0], [z]]), z, C_)

MuPAD 能否以一组向量或类似方式给出解决方案？为了使用上面的答案，我需要从那组解决方案中整理出解决方案。有没有聪明的方法来做到这一点？到目前为止，我的解决方案是找到我知道会出现在解决方案中的标志，例如'([['，然后选择下面的数字，这真的很难看，如果解决方案由于某种原因看起来与我的情况有点不同''ve涵盖了它不起作用。

编辑

当我使用@horchler 在下面的答案中建议的解决方案时，我得到的解决方案与我之前的实现相同。但是对于某些情况（不是全部），它需要更长的时间。例如。对于下面的 T，下面建议的解决方案需要一分钟以上，而使用 evalin（我以前的实现）需要一秒钟。

T = [2.4336 1.4309 0.5471 0.0934 9.5838 -0.1013 -0.2573 2.4830 ...
     36.5464 0.4898 -0.5383 61.5723 1.7637 36.0816 11.8262]

新功能：

function ut1 = testMupadSolver(T)
% # Input T should be a vector of 15 elements
allSolutions = feval(symengine,'numeric::polysysroots', ...
    [eq1(T),eq2(T)],'[u,v]');
end
function eq = eq1(T)
syms u v
eq = -T(13)*u^4 + (4*T(15) - 2*T(10) - T(11)*v)*u^3 + (3*T(13) - 3*T(6) ...
    + v*(3*T(14) -2*T(7)) - T(8)*v^2)*u^2 + (2*T(10) - 4*T(1) + v*(2*T(11) ...
    - 3*T(2)) + v^2*(2*T(12) - 2*T(3)) - T(4)*v^3)*u + v*(T(7) + T(8)*v ...
    + T(9)*v^2) + T(6);
end


function eq = eq2(T)
syms u v
eq = (T(14) - T(13)*v)*u^3 + u^2*(T(11) + (2*T(12) - 2*T(10))*v ...
    - T(11)*v^2) + u*(T(7) + v*(2*T(8) - 3*T(6) ) + v^2*(3*T(9) - 2*T(7)) ...
    - T(8)*v^3) + v*(2*T(3) - 4*T(1) + v*(3*T(4) - 3*T(2)) + v^2*(4*T(5) ...
    - 2*T(3)) - T(4)*v^3) + T(2);
end

是否有充分的理由说明为什么需要这么长时间？

Answer 1

首先，Matlab 通过字符串命令与 MuPAD 通信，因此最终无法绕过字符串的使用。而且因为它是本机格式，如果您将大量数据传递到 MuPAD，最好的方法是将所有内容快速有效地转换为字符串（sprintf 通常是最好的）。但是，就您而言，我认为您可以使用feval 而不是evalin，它允许您传入常规的Matlab 数据类型（在后台sym/feval 进行字符串转换并调用evalin）。此方法在此MathWorks article 中进行了讨论。可以使用以下代码：

T = [0 0 0 0 0 0 0 0 0 0 1 0 1 0 1];
syms u v;
eq1 = -T(13)*u^4 + (4*T(15) - 2*T(10) - T(11)*v)*u^3 + (3*T(13) - 3*T(6) ...
    + v*(3*T(14) -2*T(7)) - T(8)*v^2)*u^2 + (2*T(10) - 4*T(1) + v*(2*T(11) ...
    - 3*T(2)) + v^2*(2*T(12) - 2*T(3)) - T(4)*v^3)*u + v*(T(7) + T(8)*v ...
    + T(9)*v^2) + T(6);
eq2 = (T(14) - T(13)*v)*u^3 + u^2*(T(11) + (2*T(12) - 2*T(10))*v ...
    - T(11)*v^2) + u*(T(7) + v*(2*T(8) - 3*T(6) ) + v^2*(3*T(9) - 2*T(7)) ...
    - T(8)*v^3) + v*(2*T(3) - 4*T(1) + v*(3*T(4) - 3*T(2)) + v^2*(4*T(5) ...
    - 2*T(3)) - T(4)*v^3) + T(2);
allSolutions = feval(symengine, 'numeric::polysysroots',[eq1,eq2],'[u,v]');

最后一个参数仍然需要是一个字符串（或省略），并且将==0 添加到等式中也不起作用，但无论如何，零是隐含的。

对于第二个问题，numeric::polysysroots 返回的结果很不方便，不好处理。这是一组 (DOM_SET) 矩阵。我尝试使用coerce 将结果转换为其他内容，但无济于事。我认为您最好将输出转换为字符串（使用char）并解析结果。我这样做是为了更简单的输出格式。我不确定它是否会有所帮助，但请随意查看我的 sym2float，它仅使用一些优化处理符号矩阵（'matrix([[ ... ]])' 部分是您的输出）。

最后一件事。您的辅助函数是否包含多余的括号？这似乎足够了

sT = @(x)num2str(T(x),17);

或

sT = @(x)sprintf('%.17g',T(x));

注意num2str 默认只转换为四位小数。 int2str（如果T(x) 始终为整数，则应使用%d）。