empty sym:0-by-1 error nonlinear equations
2 次查看(过去 30 天)
显示 更早的评论
Hello,
i have 5 nonlinear equations and i wrote code;
clear all;
clc;
syms u1 u2 u3 u4 u5;
eq1=u1== 0.5 *(0.2*u2)+(u1^3);
eq2=0.2*u2==0.8*u1+0.2*u3 +(u2^3);
eq3=0.2*u3== 0.5*(0.8*u2+0.2*u4)+(u3^3);
eq4=0.2*u4==0.8*u3+0.2*u5 +(u4^3);
eq5=0.2*u5== 0.5*(0.8*u4)+(u5^3);
eqs=[eq1,eq2,eq3,eq4,eq5];
[u1,u2,u3,u4,u5]=vpasolve(eqs,[u1,u2,u3,u4,u5])
it works.
But if i have 6 nonlinear equation and i wrote code,
clear all;
clc;
syms u1 u2 u3 u4 u5 u6;
eq1=0.2*u1== 0.5 *(0.2*u2)+(u1^3);
eq2=0.2*u2==0.8*u1+0.2*u3 +(u2^3);
eq3=0.2*u3== 0.5*(0.8*u2+0.2*u4)+(u3^3);
eq4=0.2*u4==0.8*u3+0.2*u5 +(u4^3);
eq5=0.2*u5== 0.5*(0.8*u4+0.2*u6)+(u5^3);
eq6=0.2*u6==0.8*u5 +(u6^3);
eqs=[eq1,eq2,eq3,eq4,eq5,eq6];
[u1,u2,u3,u4,u5,u6]=vpasolve(eqs,[u1,u2,u3,u4,u5,u6])
it does not work and i get a error like empty sym:0-by-1. Could you please help me?
Do you have any suggestions make it easier?
Thanks in advance.
2 个评论
Walter Roberson
2020-12-19
The first 5 equations have solutions as a system of 243 values involving polynomials of degree 243 in u6.
Those solutions can be computed (in placeholder form as individual roots) and substituted into the 6th equation, and then try to find solutions. If there are solutions, you would expect 729 of them (possibly not unique). It is likely that the majority of them would be complex valued.
回答(2 个)
Kiran Felix Robert
2020-12-18
编辑:Kiran Felix Robert
2020-12-18
Hi Semiha,
vpasolve returns zero when the system of equations does not have a solution, refer tips section of vpasolve documentation.
0 个评论
Walter Roberson
2020-12-21
u = sym('u', [1 6]); syms(u);
eqn = [1/5*u1 == 1/5*1/2*u2 + u1^3, 1/5*u2 == 4/5*u1 + 1/5*u3 + u2^3, 1/5*u3 == 1/2*(4/5*u2 + 1/5*u4) + u3^3, 1/5*u4 == 4/5*u3 + 1/5*u5 + u4^3, 1/5*u5 == 1/2*(4/5*u4 + 1/5*u6) + u5^3, 1/5*u6 == 4/5*u5 + u6^3]
sol = solve(eqn(1:5), u(1:5));
E6 = subs(eqn(6), sol);
R6 = findSymType(E6, 'root');
ch6 = children(R6);
co6 = coeffs(ch6{1}, ch6{2}, 'all');
CO6 = matlabFunction(co6);
syms R
E6C = subs(lhs(E6)-rhs(E6), R6, R)
E6CF = matlabFunction(E6C);
FUNraw = @(u6) E6CF(roots(CO6(u6)), u6);
norm_sq = @(x) x.*conj(x);
FUN = @(u6) norm_sq(FUNraw(u6))
U6 = linspace(-10,10);
values = arrayfun(FUN, U6, 'uniform', 0);
vmat = cell2mat(values);
[bestval, bestidx] = min(vmat,[],'all', 'linear');
[row,col] = ind2sub(size(vmat), bestidx);
fprintf('Best root is #%d\n', row)
fprintf('Best residue found is %g near u6 = %g\n', bestval, U6(col))
fprintf('\nRefining...\n\n');
refU6 = linspace(U6(col-1), U6(col+1), 1000);
refvalues = arrayfun(FUN, refU6, 'uniform', 0);
refvmat = cell2mat(refvalues);
[refbestval, refbestidx] = min(refvmat,[],'all', 'linear');
[refrow,refcol] = ind2sub(size(refvmat), refbestidx);
bestrefU6 = refU6(refcol);
fprintf('Best refined root is #%d\n', refrow)
fprintf('Best refined residue found is %g near u6 = %g, in [%g..%g]\n', refbestval, bestrefU6, refU6(refcol-1), refU6(refcol+1))
plot(refU6(refcol-10:refcol+10), refvmat(refrow, refcol-10:refcol+10))
goodu1 = double(subs(sol.u1, u6, bestrefU6));
goodu2 = double(subs(sol.u2, u6, bestrefU6));
goodu3 = double(subs(sol.u3, u6, bestrefU6));
goodu4 = double(subs(sol.u4, u6, bestrefU6));
goodu5 = double(subs(sol.u5, u6, bestrefU6));
goodu6 = double(bestrefU6);
fprintf('u1 might be near %g %+gI\n', real(goodu1), imag(goodu1) );
fprintf('u2 might be near %g %+gI\n', real(goodu2), imag(goodu2) );
fprintf('u3 might be near %g %+gI\n', real(goodu3), imag(goodu3) );
fprintf('u4 might be near %g %+gI\n', real(goodu4), imag(goodu4) );
fprintf('u5 might be near %g %+gI\n', real(goodu5), imag(goodu5) );
fprintf('u6 might be near %g %+gI\n', real(goodu6), imag(goodu6) );
sol = vpasolve(eqn, u, [goodu1; goodu2; goodu3; goodu4; goodu5; goodu6])
if isempty(sol)
fprintf('... but vpasolve cannot find a solution near there :(\n');
end
u1 might be near -0.011815 +0I
u2 might be near -0.0236135 +0I
u3 might be near 0.0237123 +0I
u4 might be near 0.141745 +0I
u5 might be near 0.0326565 +0I
u6 might be near -0.502017 +0I
This code numerically finds a position that the equations almost hold, then refines it once in order to find a better precision. It reports on what it finds, and then tries to use vpasolve() starting from there... but will probably fail the vpasolve.
With the complexity of the system, it is difficult to say for sure whether any solution for all of the variables exists. The first five equations can be solved for the first five variables, giving you a polynomial of degree 243 in u6. The code then substitutes in the polynomial into the 6th equation, and then tries to minimize the absolute value difference between the left and right side of the 6th equation.
The code does not assume that particular roots out of the 243 are going to be the only valid roots: for each potential u6 value, the code evaluates all 243 roots numerically and uses the results in the 6th equation. The "best" value is determined by the norm-squared of the complex residues, norm(A+B*1i)^2 = A^2 + B^2 which just happens to be calculated by (A+B*1i).*conj(A+B*1i)
All of this is done numerically because matlabFunction() is not able to generate code for the internal RootOf that are generated for the 243 roots of the first five equations.
The code does, though, assume that u6 is real valued and somewhere in the range -10 to +10, which was a guess on my part. The location it finds as a decently low residue, so the guess was not out of line. It could be the case that there is an actual solution outside that range, possibly requiring a complex u6. Or it might be the case that there is no solution at all. We reduced it down to a single complicated equation involving roots of a degree 243 polynomial; if we could solve the last equation then the remaining variables would fall out of that known u6.
4 个评论
Walter Roberson
2021-1-2
With even just 6 equations, the numbers in the polynomials get large enough that you cannot work in floating point without too much loss of precision. I think you would have a quite difficult time doing the calculations.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Polynomials 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!