Solving system non-linear equations

Dear Matlab Community,

I would like to ask a question regarding the solving of a 4 equations non-linear system.

I have tried two approaches: the first one by defining a system and use solve to extract the solutions, but I get into an infinite loop and Matlab doesn't give any solution:

syms E1 X P E2 
eq1 = E0 + E1*(X)^(sr(1))+E2*(P)^(sr(1)) == Eapp(1);
eq2 = E0 + E1*(X)^(sr(2))+E2*(P)^(sr(2)) == Eapp(2);
eq3 = E0 + E1*(X)^(sr(3))+E2*(P)^(sr(3)) == Eapp(3);
eq4 = E0 + E1 + E2 == Ei
eqns = [eq1 eq2 eq3 eq4];
vars = [E1 X P E2];
A = solve(eq1, eq2, eq3, eq4)

My unknown are E1, X, E2, P and all the others are defined parameters.

I have also tried to implement a function:

function F = root2d(x)

sr = [0.05 0.1 0.5 1];
Eapp = [2.09 2.18 2.27 2.4].*1e6; 
Ei = 2.55e6;
E0 = 2.14e6;

F(1) = E0 + x(1)*(x(2))^(sr(1))+x(3)*(x(4))^(sr(1)) - Eapp(1);
F(2) = E0 + x(1)*(x(2))^(sr(2))+x(3)*(x(4))^(sr(2)) -Eapp(2);
F(3) = E0 + x(1)*(x(2))^(sr(3))+x(3)*(x(4))^(sr(3)) - Eapp(3);
F(4) = E0 + x(1) + x(3) - Ei;

end

Said function is called in the main script and I have used fsolve to extract the solutions:

fun = @root2d;
x0 = [1.5e6,1,1.5e6,1];
x = fsolve(fun,x0)

In this case, my problem is that I cannot properly set an initial interval x0 where I can look for the solutions, since I cannot plot a 4 variables function to visually see the zeros, I get the following error (I have estimated the values of E1 and E2 as order of magnitude but I have no clue of which the values of P and X should be):

No solution found.

fsolve stopped because the relative size of the current step is less than the

value of the step size tolerance squared, but the vector of function values

is not near zero as measured by the value of the function tolerance.

Does anyone have any tip to solve this problem?

Lucrezia

[Walter Roberson]

My tests suggest that there is no solution. You can get formulas in terms of roots of degree 12 equations, but when you back-substitute through, the results are inconsistent.

It is the sort of situation where you can solve one side of an equation to -1 and the other side to sqrt(-1) and your processing had thought that you had found a complete solution because your processing had involved taking the 4th power of both side.

[Lucrezia]

Hi Walter,

Thank you for your reply.

It is quite surprisingly since these are textbook and known equations. I will have a closer look if there is any inconsistence related to order of magnitude or numerical values.

Lucrezia

Answer 1

[Matt J]

You could do a numerical grid search using ndgrid to see where if anywhere, the approximate roots lie within a broad set of bounds. If there are any roots, you can then do a more refined search with lsqnonlin in the neighborhood of those.

In the code below, I've sketched this, although note that I have used the 4th equation to eliminate x1 and thereby reduce the amount of grid data that is needed.

[X2,X3,X4]=ndgrid( linspace(a,b,500), linspace(c,d,500) , linspace(e,f,500) );

F=gridsearch(X2,X3,X4);

[minval,location]=min(F(:))


function F = gridsearch(x2,x3,x4)

    sr = [0.05 0.1 0.5 1];
    Eapp = [2.09 2.18 2.27 2.4].*1e6; 
    Ei = 2.55e6;
    E0 = 2.14e6;

    F1 = E0 + (Ei-E0-x3).*(x2).^(sr(1))+x3.*(x4).^(sr(1)) - Eapp(1);
    F2 = E0 + (Ei-E0-x3).*x2.^(sr(2))+x3.*(x4).^(sr(2)) - Eapp(2);
    F3 = E0 + (Ei-E0-x3).*(x2).^(sr(3))+x3.*(x4).^(sr(3)) - Eapp(3);

    F=abs(F1)+abs(F2)+abs(F3);

end

[Lucrezia]

Hi Matt,

Thank you very much for the useful suggestion. I have run the script and set initial values:

a=0, b=500, c=0, d=500, e=0, f=500.

(First question: how should I set a,b,c,d,e,f "a priori"?)

As far as I understood, said variables define a 3 dimensional space (for F1, F2, F3) where I can look for my roots for F1, F2 and F3, respectively.

The value Matlab returns me with these values is:

Location: 124999501

minVal: 2.085710479943443e+05

I have some questions regarding this result:

1. Should the location be in the interval I have defined within linspace?
2. I would have expect also a different root for x2 and x4 (as order of magnitude), but I got only one minimum: how should I use this result in my original script, in the part where I use fsolve and the zeros? Even if I reduce the system to three equations, I would still need to use as inputs three zeros for x0 (fsolve(fun, x0)).

Thank you for the clarifications,

Lucrezia

[Matt J]

(First question: how should I set a,b,c,d,e,f "a priori"?)

Based on your knowledge of the physics of the problem, you should choose the grid bounds to cover a region the root(s) are likely to lie within. Note that with,

a=0, b=500, c=0, d=500, e=0, f=500.

you are saying you believe solutions will be found in the interval 0<= x2,x3,x4 <=500. In the code you originally posted, you were originally guessing a solution of x0 = [1.5e6,1,1.5e6,1] which is well outside this interval.

Should the location be in the interval I have defined within linspace?

Yes. to obtain the x2,x3,x4 coordinates where the minimum was found, you would do,

x2=X2(location);
x3=X3(location);
x4=X4(location);

I would have expect also a different root for x2 and x4 (as order of magnitude), but I got only one minimum:

The min() command only reports the first global minimum it encounters. The purpose of what I presented was to see, as a first step, if minval is sufficiently close to zero to believe that there is a root within the grid box that you have chosen. To me, it doesn't look very close to zero, and Walter has suggested that you won't be able to find one for any box that you choose. Nevertheless, if you have faith that this minval represents a root, you can collect additional approximate roots by doing,

region=(F<=minval+tolerance);

x2=X2(region);
x3=X3(region);
x4=X4(region);

Answer 2

[Walter Roberson]

syms x [1 4]

sr = sym([0.05 0.1 0.5 1]);
Eapp = sym([2.09 2.18 2.27 2.4]).*1e6; 
Ei = sym(2.55).*1e6;
E0 = sym(2.14).*1e6;

F(1) = E0 + x(1)*(x(2))^(sr(1))+x(3)*(x(4))^(sr(1)) - Eapp(1);
F(2) = E0 + x(1)*(x(2))^(sr(2))+x(3)*(x(4))^(sr(2)) -Eapp(2);
F(3) = E0 + x(1)*(x(2))^(sr(3))+x(3)*(x(4))^(sr(3)) - Eapp(3);
F(4) = E0 + x(1) + x(3) - Ei;

obj = simplify(expand(sum(F.^2)))

obj = 

crit_x1 = solve(diff(obj,x(1)),x(1))

crit_x1 = 

obj2 = subs(obj, x(1), crit_x1)

obj2 = 

crit_x3 = solve(diff(obj2,x(3)),x(3))

crit_x3 = 

obj3 = subs(obj2, x(3), crit_x3)

obj3 = 

Unfortunately, MATLAB is not able to solve obj3 for x(2) or x(4).

Now let us look more closely at obj3

[N,D] = numden(obj3);

D

D = 

With the fractional powers, we can rule out negative x(2) and x(4) for all practical purposes. So each individual term is positive and there are no subtractions, and there is a constant +1 at the end so the value cannot be 0. Therefore the denominator is non-zero for all values we care about (though there are potentially complex x(2), x(4) that would lead to a denominator of 0). Non-zero denominator means we do not need to worry about singularities caused by the denominator, and we can concentrate on the numerator

N

N = 

That does have subtractions, so we have the possibility of zeros.

fsurf(N, [0 2 0 2]); view(3)

So, possibly along 0

x2lim = limit(N, x(2), 0, 'right')

x2lim = 

x4lim = limit(N, x(4), 0, 'right')

x4lim = 

figure(); fplot(x2lim, [0 2])

figure(); fplot(x4lim, [0 2])

limit(x2lim, x(4), 0, 'right')

ans = 

18500000000

limit(x4lim, x(2), 0, 'right')

ans = 

18500000000

So... unless the shape changes quite a bit with higher values, or the calculas solutions are for maxima instead of minima, or you go to complex calculations, there is no solution to the set of equations.