3

我的问题是是否有在 Matlab 脚本中使用 MuPAD 函数的好方法。背景是我有一个问题，我需要找到一组非线性方程的所有解。以前的解决方案是solve在 Matlab 中使用，它适用于我的一些模拟（即一些输入集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 可以将解决方案作为一组向量或类似方式给出吗？为了使用上面的答案，我需要从那组解决方案中整理出解决方案。有没有聪明的方法来做到这一点？到目前为止，我的解决方案是找到我知道将出现在解决方案中的迹象，例如'([['并选择下面的数字，这真的很难看，如果由于某种原因解决方案看起来与我所涵盖的情况有点不同它不起作用。

编辑

当我使用@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

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

4

1 回答 1

2

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

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（或者如果始终是整数%d，则应该使用）。T(x)

于 2013-12-31T07:04:58.720 回答