Поиск всех решений системы нелинейных уравнений с помощью MuPAD

Мой вопрос в том, есть ли хороший способ использовать функции MuPAD в сценарии Matlab. Предыстория заключается в том, что у меня есть проблема, в которой мне нужно найти все решения набора нелинейных уравнений. Предыдущее решение заключалось в использовании 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

Есть ли веская причина, почему это занимает так много времени?


matlab symbolic-math mupad
person Fija    schedule 30.12.2013    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Во-первых, 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);

or

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

Обратите внимание, что num2str по умолчанию преобразуется только в четыре десятичных разряда. int2str (или %d следует использовать, если T(x) всегда является целым числом).

person horchler    schedule 31.12.2013
comment
Привет и спасибо! Я начну с ответа на ваш последний комментарий: лишние скобки нужны, поскольку элементы в T не обязательно положительны, поэтому в моей текущей версии они необходимы, но, возможно, они не соответствуют вашим предложениям. Что касается других ваших предложений, они выглядят действительно хорошо, но я не могу проверить их должным образом, пока не вернусь на работу через пару дней. - person Fija; 31.12.2013
comment
пс. У меня слишком мало репутации, чтобы голосовать за ваш ответ, но я сделаю это, как только смогу. д.с. - person Fija; 31.12.2013
comment
Пожалуйста, смотрите мое редактирование @horchler. Я реализовал ваше предложение, и оно отлично работает, за исключением того, что иногда это занимает гораздо больше времени. - person Fija; 03.01.2014
comment
@Fija: Вы уверены, что evalin не займет столько же времени? Символический набор инструментов кэширует результаты, поэтому, если вы сначала запустите feval для нового уравнения, вы можете не получить справедливого сравнения. Вы можете попробовать позвонить clear all и/или clear classes. Могут быть и другие вещи. В любом случае вы можете ожидать, что feval будет немного медленнее по причинам, изложенным в моем ответе. - person horchler; 03.01.2014
comment
Нет, я уверен, что выполнял и clear all, и clear classes между прогонами, и feval занимает примерно 70 секунд, а evalin примерно 2 секунды. Но я думаю, что ваш ответ был хорошим, и это другая проблема, поэтому я приму ваш ответ. знак равно - person Fija; 04.01.2014
comment
@Fija: Может случиться так, что в некоторых случаях на промежуточном этапе могут быть возвращены сложные аналитические решения, которые дорого оценивать численно. Возможно, что то, как вы указываете код evalin, требует непосредственной оценки числового решения (символический класс, но полностью числовое). Последнее может быть намного быстрее, если формулы, возвращаемые feval, являются дорогостоящими. Вы могли бы видеть, можете ли вы использовать vpa на T или ваших уравнениях. - person horchler; 04.01.2014
comment
Я предполагаю, что это может быть и наоборот, и что в некоторых случаях оценивается числовое решение, которое оказывается медленным. Дело в том, что строка, которую создает feval, вероятно, отличается числовыми константами/параметрами. Включение/исключение десятичной точки (.0) может изменить внутреннюю оценку вещей. - person horchler; 04.01.2014