empty sym:0-by-1 error nonlinear equations

5 views (last 30 days)
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 Comments
Alex Sha
Alex Sha on 19 Dec 2020
Hi, zero for all u are the solution
Walter Roberson
Walter Roberson on 19 Dec 2020
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.

Sign in to comment.

Answers (2)

Kiran Felix Robert
Kiran Felix Robert on 18 Dec 2020
Edited: Kiran Felix Robert on 18 Dec 2020
Hi Semiha,
vpasolve returns zero when the system of equations does not have a solution, refer tips section of vpasolve documentation.
In this case you can provide an initial guess to help the solver find a solution.
For solving non-linear equation you can also try fsolve.

Walter Roberson
Walter Roberson on 21 Dec 2020
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 Comments
Semiha Tombuloglu
Semiha Tombuloglu on 2 Jan 2021
What about if i want to find u20 ?
I tried both Mathematica and Matlab to get solution. It is always running.
I don't know if it has an easier solution.
Walter Roberson
Walter Roberson on 2 Jan 2021
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!