empty sym:0-by-1 error nonlinear equations

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.

[Alex Sha]

Hi, zero for all u are the solution

[Walter Roberson]

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.

Answer 1

[Kiran Felix Robert]

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.

Answer 2

[Walter Roberson]

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.

[Semiha Tombuloglu]

We can suppose that all the unknowns (u1,u2, u3, u4, u5, u6) are real.

Actually i need more equations and more unknowns. Do you think that finding these values depend on computer speed? Is it really imposible to find?

[Walter Roberson]

If you solve the first 5 equations, you get values in terms of the roots of a polynomial of degree 81. If you replace that root with a variable, you get a 6th equation that can be solved in terms of that variable without taking any more roots; powers of that temporary variable up to 161 are involved. You can then ask when the result would be real-valued. It turns out that there is a range of values for the root between approximately -2162/10000 and 299063/500000 where it would be real-valued -- but not that whole range, with a section approximately 0.025 to 0.225 where the expression for u6 would be imaginary valued.

For example the expression is real-values when the root is 1/2 .

You can take any value in that range and substitute it in as the root of the polynomial of degree 81 to get an expression in terms of u6 and that expression would need to be 0

RootOf(9313225746154785156250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*z^243 - 150874257087707519531250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*z^241 and so on + u6, z)

substitute in the known z value (a value in the range -2162/10000 and 299063/500000 excluding the center part mentioned) to solve for u6. For example when the root is 1/2 then

u6 = 6580823775028600461390155895364563919719373134399899714673415325455 / 5708990770823839524233143877797980545530986496

which is about 1.15E21 .

Substitute that and the picked root back into the solution to the first five equations in order to get solutions for the whole system.

Well... in theory. In practice I am getting contradictions when I back substitute.:(

[Semiha Tombuloglu]

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]

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.