vpasolve Больше неизвестных, чем уравнений => Нет решений - PullRequest
Новая
       13

vpasolve Больше неизвестных, чем уравнений => Нет решений

0 голосов
/ 21 мая 2019

У меня есть 3 уравнения и 4 неизвестных, которые я хочу решить. Конечно, есть решения, но vpasolve не возвращает никаких, если я опускаюсь до 3 экв. И 3 неизвестных, это работает хорошо. Я знаю, что с большим количеством неизвестных у меня есть почти бесконечное количество решений, так как мне заставить это работать в этом случае?

Вот мой код: 

syms x y z theta1 theta2 theta3 phi1

xEquation = x == cos(theta1)*cos(phi1) + cos(theta1 + theta2)*cos(phi1) + cos(theta1 + theta2 + theta3)*cos(phi1)
yEquation = y == cos(theta1)*sin(phi1) + cos(theta1 + theta2)*sin(phi1) + cos(theta1 + theta2 + theta3)*sin(phi1)
zEquation = z == sin(theta1) + sin(theta1 + theta2) + sin(theta1 + theta2 + theta3)

x = 2;
y = 1.5;
z = 0;
sol = vpasolve([eval(xEquation), eval(yEquation), eval(zEquation)], [theta1, theta2, theta3, phi1], [-pi pi; -pi pi; -pi pi; -pi pi;]);

Это производит структуру sol с 4 полями, но они пусты, без решений.

1 Ответ

1 голос
/ 21 мая 2019

Решение m уравнений с n неизвестными, такими как m<n, означает, что некоторые переменные будут параметрами, зависящими от других переменных.

Например 

x-3y+z = 2
x-y+5z = 5

Предположим, z является параметром
Чтобы решить эту проблему в Matlab вот код 

syms x y z

eq1 = x-3*y+z ==2;
eq2 = x-y+5*z ==5;
sol = solve(eq1, eq2, x, y);
sol.x
sol.y

Как видите, z опущено в выражении решения, что означает, что оно будет рассматриваться как параметр
Решение 

sol.x = 13/2 - 7*z
sol.y = 3/2 - 2*z
  • Из решения видно, что x, y and z не являются числовыми значения, поэтому вы не можете использовать vpasolve, который обозначает Variable-precision arithmetic. vpasolve дает точность численного решения
  • Также, поскольку x, y and z не являются числовыми значениями, вы не можете предопределить диапазон для них, если вы не исправите z first

Вы можете использовать solve для просмотра набора решений, здесь я рассматриваю phi1 для параметра, поэтому опускаю его в выражении решения 

syms x y z theta1 theta2 theta3 phi1

xEquation = 2 == cos(theta1)*cos(phi1) + cos(theta1 + theta2)*cos(phi1) + cos(theta1 + theta2 + theta3)*cos(phi1);
yEquation = 1.5 == cos(theta1)*sin(phi1) + cos(theta1 + theta2)*sin(phi1) + cos(theta1 + theta2 + theta3)*sin(phi1);
zEquation = 0 == sin(theta1) + sin(theta1 + theta2) + sin(theta1 + theta2 + theta3);


sol = solve(xEquation, yEquation, zEquation, theta1, theta2, theta3);

Набор решений 

sol.theta1 = [-pi*(2*n - 1); pi*(2*m + 1); pi*k; pi*k]
sol.theta2 = [pi*(2*m + 1);  -pi*(2*n - 1);  -pi*(2*n - 1); z]
sol.theta3 = [pi*k; pi*k; pi*(2*m + 1); pi*(2*m + 1)] 
phi1 is the parameter

Первый набор решений 

X = [sol.theta1(1); sol.theta2(1); sol.theta3(1); phi1]
X = [-pi*(2*n - 1); pi*(2*m + 1); pi*k; phi1]

В соответствии с вышеуказанными вычетами написано 4 комплекта

  • z - параметр, k, m, n - целые числа, которые в основном используется для периодичности тригонометрических функций

  • Если вы установите z в диапазоне [-pi, pi], вы можете настроить k, m and n на получить правильное решение в диапазоне [-pi, pi].

Как видите, время стоит

Использование fmincon

  • В качестве альтернативы вы можете использовать fmincon для решения вашей проблемы
  • Я в основном минимизирую постоянную функцию, fmincon просто выдаст решения, которые удовлетворяют ограничениям, здесь ceq = 0
  • Интервал поиска [-pi pi] преобразуется в lb = -pi и ub = pi
  • Я установил начальное предположение на 0

Код следующий 

t = 0:0.1:1;
x = 1.5 + 0.5 .* cos(8 .* pi .* t);
y = 1.5 + 0.5 .* sin(8 .* pi .* t); 
z = 1 .* t .* ones(size(x));

lb = -pi*ones(1, 4);
ub = -lb;

p0 = zeros(1,4);
sol = cell(1,length(t));

for i = 1:length(t)
    sol{i} = fmincon(@(p)0,p0,[],[],[],[],lb,ub,@(p)nonlincon(x(i),y(i), z(i), p(1), p(2), p(3), p(4)));


end


function [c, ceq] = nonlincon(x,y, z, theta1, theta2, theta3, phi1)

    c = [];
    ceq(1) = cos(theta1)*cos(phi1) + cos(theta1 + theta2)*cos(phi1) + cos(theta1 + theta2 + theta3)*cos(phi1)-x;
    ceq(2) = cos(theta1)*sin(phi1) + cos(theta1 + theta2)*sin(phi1) + cos(theta1 + theta2 + theta3)*sin(phi1)-y;
    ceq(3) = sin(theta1) + sin(theta1 + theta2) + sin(theta1 + theta2 + theta3)-z;
end

Решение

Второй набор решений, когда t = 0.1 равен sol{2}

sol{2}.(1) = pheta1
sol{2}.(2) = pheta2
sol{2}.(3) = pheta3
sol{2}.(4) = phi1

Вы можете следовать той же логике, чтобы найти решение в разное время t

Целое решение

enter image description here

...