MATLAB Answers

Solving system non-linear equations

18 views (last 30 days)
Lucrezia Maini
Lucrezia Maini on 13 May 2021
Answered: Walter Roberson on 14 May 2021
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
  2 Comments
Lucrezia Maini
Lucrezia Maini on 14 May 2021
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

Sign in to comment.

Accepted Answer

Matt J
Matt J on 13 May 2021
Edited: Matt J on 13 May 2021
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
  2 Comments
Matt J
Matt J on 14 May 2021
(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);

Sign in to comment.

More Answers (1)

Walter Roberson
Walter Roberson on 14 May 2021
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.

Community Treasure Hunt

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

Start Hunting!