solving a non-linear problem
Show older comments
Hello, Below is a problem I need to solve:
System (1) of 7 equations (1a,b,c,d,e,f), 7 unknowns x(1...7), and 4 parameters a,n,k and w;
1a. x(1)/x(5)=n
1b. x(2)/x(5)=k
1c. x(4)^2*x(3)*x(5)=a
1d. x(4)^2*x(6)*x(5)=w
1e. x(1)+x(2)+x(5)-x(3)-x(6)=0
1f. x(7)=1/2*(x(1)+x(2)+x(3)+x(5)+x(6))
1g. x(4)=10^(-x(7)^2/(1+x(7)^2))
with a=10^(-9.9727); k=10^(3.8999); n=10^(4.6711); w=8.472e-11 xi must be real and positive, initial guess may be x(1)=0.1, x(2)=0.001, x(3)=0.001, x(4)=1, x(5)=10e-7, x(6)=0.1, x(7)=0.1
Ultimately, I need a script to solve 1000 such system with Xi are matrices of unknowns and a,k,n and w are matrices of parameters. So I would be very grateful if someone could suggest to me a convenient script to do it systematically. I have tried with fsolve so far but cannot reach a satisfactory convergence for some reasons.
Thank you in advance.
Accepted Answer
Roger Stafford
on 17 Oct 2014
The most recent version of your eight equations (9/11/14) in which eq. 6. had been changed was this:
1. x1/x5 = n;
2. x6*x2/(x7*x5)^2 = c
3. x7^2*x3*x5 = w
4. x6*x4*(x7*x5)^2 = co
5. x1+2*x2+x5-x3-2*x4 = 0
6. 1/2*(x1+4*x2+x5+x3+4*x4) = x8
7. x7 = 10^(-x8^(1/2)/(1+x8^(1/2))+0.2*x8)
8. x6 = 10^(-4*x8^(1/2)/(1+x8^(1/2))+0.8*x8)
Looking ahead to avoid a quartic equation, if eqs. 5. & 6. hold, they may be combined to form the equation
9. 1/2*x1+1/2*x5-3/2*x3-4*x4+x8 = 0
If a value of x8 is given, then x7 and x6 can be obtained from 7. & 8. Then from 1., 3., and 4. it follows that the cubic equation in x5 must hold
(n+1)*x5^3+2*x8*x5^2-3*w/x7^2*x5-8*co/x6/x7^2 = 0
where x5 is the only unknown. Then x5 can be determined by
r = roots([n+1,2*x8,-3*w/x7^2,-8*co/x6/x7^2]);
x5 = r(imag(r)==0 & real(r)>0);
(since it can be shown that because of the signs in the above cubic it must have exactly one real positive root.) All other x's can then be easily computed. Therefore we can use 'fzero' to adjust values of x8 so as to obtain a zero value in eq. 5.
The following is code that should accomplish this. Let n, c, w, and co be defined. Write the subfunction:
function z = myfun(x8)
t = x8.^(1/2);
x7 = 10.^(-t./(1+t)+0.2*x8);
x6 = 10.^(-4*t./(1+t)+0.8*x8);
x5 = zeros(size(x8));
for k = 1:length(x8)
r = roots([n+1,2*x8(k),-3*w/x7(k)^2,-8*co/x6(k)/x7(k)^2]);
x5(k) = r(imag(r)==0 & real(r)>0);
end
x1 = n*x5;
x2 = c./x6.*(x7.*x5).^2;
x3 = w./x7.^2./x5;
x4 = co./x6./(x7.*x5).^2;
z = x1+2*x2+x5-x3-2*x4;
return
and call on it as follows:
x8 = fzero(@myfun,x8est);
t = x8^(1/2);
x7 = 10^(-t/(1+t)+0.2*x8);
x6 = 10^(-4*t/(1+t)+0.8*x8);
r = roots([n+1,2*x8,-3*w/x7^2,-8*co/x6/x7^2]);
x5 = r(imag(r)==0 & real(r)>0);
x1 = n*x5;
x2 = c/x6*(x7*x5)^2;
x3 = w/x7^2/x5;
x4 = co/x6/(x7*x5)^2;
where 'x8est' has been carefully chosen as an initial estimate for which successful convergence by 'fzero' can be achieved. (Not every estimate value will succeed.)
More Answers (6)
Roger Stafford
on 16 Aug 2014
Edited: Roger Stafford
on 16 Aug 2014
3 votes
In spite of your claim of seven equations you have given only six of them, which are not enough to uniquely determine the seven unknowns. If that is what you gave to 'fsolve', it would not be happy. There would be infinitely many solutions. Even with the constraint that all x(i) be positive it is likely there are still infinitely many solutions.
Also you have not used the parameter w in any of these six equations. It leads me to suspect you have a seventh equation which you have neglected to mention.
Roger Stafford
on 17 Aug 2014
By appropriate algebraic manipulations you can reduce your (corrected) seven equations to finding the solution of a single equation in x5:
10^(-((n+k+1)*x5)^2/(1+((n+k+1)*x5)^2)) - sqrt((a+w)/(n+k+1))/x5 = 0
The other x values can then be found in terms of x5:
x4 = sqrt((a+w)/(n+k+1))/x5;
x1 = n*x5;
x2 = k*x5;
x3 = a/x5/x4^2;
x6 = w/x5/x4^2;
x7 = 1/2*(x1+x2+x3+x5+x6);
If the four parameters, a, w, n, and k are all positive, it follows that all x values will be positive.
Therefore you can use matlab's fzero function which is less complicated than calling on fsolve. You will need to experiment to find the best estimate x0 for solving for x5 with fzero. It just has to be in the right approximate neighborhood and may depend to some extend on the values of the four parameters a, w, n, and k.
3 Comments
Roger Stafford
on 19 Aug 2014
The answer to 3. is very similar to your original seven equations. The reasoning goes as follows (use x1, x2, etc. instead of x(1), x(2), etc.)
x1 = n*x5
x2 = k*x5
x3 = a/x4^2/x5
x6 = w/x4^2/x5
x8 = s/x4^2/x5
x1+x2+x5 = (n+k+1)*x5
x3+x6+x8 = (a+w+s)/x4^2/x5
0 = (x1+x2+x5)-(x3+x6+x8) = (n+k+1)*x5-(a+w+s)/x4^2/x5
(x4*x5)^2 = (a+w+s)/(n+k+1)
x4 = sqrt((a+w+s)/(n+k+1))/x5
x7 = 1/2*((x1+x2+x5)+(x3+x6+x8)) = 1/2*((n+k+1)*x5+(a+w+s)/x4^2/x5)
= 1/2*((n+k+1)*x5+(a+w+s)/(x4*x5)^2*x5)
= 1/2*((n+k+1)*x5+(a+w+s)/((a+w+s)/(n+k+1))*x5)
= (n+k+1)*x5
x4 = 10^(-x(7)^(1/2)/(1+x(7)^(1/2)))
sqrt((a+w+s)/(n+k+1))/x5 = 10^(-x(7)^(1/2)/(1+x(7)^(1/2))) =
= 10^(-((n+k+1)*x5)^(1/2)/(1+((n+k+1)*x5)^(1/2)))
sqrt((a+w+s)*(n+k+1)) =
= ((n+k+1)*x5)*10^(-((n+k+1)*x5)^(1/2)/(1+ ((n+k+1)*x5)^(1/2)))
Let t = (n+k+1)*x5 and p = sqrt((a+w+s)*(n+k+1))
f(t) = t*10^(-sqrt(t)/(1+sqrt(t)))-p = 0
Use fzero to solve for t for the given value of p
f = @(t) t.*10.^(-sqrt(t)./(1+sqrt(t)))-p;
tt = fzero(f,[0,10*p]); % The solution must lie between 0 and 10*p
x5 = tt/(n+k+1);
The other x's can be evaluated in terms of x5 by the above equations.
The answer to 1. is that you cannot do your solutions with fzero in terms of matrices. It will have to be done in some kind of for-loop, one for each different value of p. (It is a fallacy to think that your computer can evaluate terms simultaneously . You will have to await quantum computers for that.)
The answer to 2. is that I don't know of any matlab function that can accomplish the above reduction to a single variable automatically. It takes a little sweat and work with a pen and paper. Fortunately there is still something we mere mortals are needed for that matlab can't do for us.
Roger Stafford
on 29 Aug 2014
With some mental effort I was able to reduce your new set of equations to a single variable, but this variable is x7, not x4. You can write a subfunction which you hand to 'fzero' as follows.
function F = Mattsfunc(x7)
t = sqrt(x7);
x5 = 10^(-t/(1+t)+0.2*x7);
x6 = x5^4;;
x4 = (-(n+1)*x5+sqrt((n+1)^2*x5^2+12*c*x6^2*x7*x5))/(6*c*x6^2);
x1 = n*x4;
x2 = c*x6^2*x4^2/x5;
x3 = x1+2*x2+x4;
F = x5^2*x3*x4-w;
return
The equation for x4 is derived as follows. (The equations or reasonings used at each step are indicated in the parentheses.)
x7 = 1/2*(x1+4*x2+x3+x4) (1e.)
= 1/2*(x1+4*x2+x1+2*x2+x4+x4) (1d.)
= 1/2*(2*(n+1)*x4+6*x2) (1a.)
= (n+1)*x4+3*c*x6^2*x4^2/x5 (1b.)
3*c*x6^2*x4^2+(n+1)*x5*x4-x7*x5 = 0 (Multiply by x5)
x4 = (-(n+1)*x5+sqrt((n+1)^2*x5^2+12*c*x6^2*x7*x5))/...
(6*c*x6^2) (<-- Solution to quadratic equation in x4)
The quadratic equation has two roots, but the one with a plus sign above needs to be chosen if x4 is to be positive.
The other steps are evident from your original equations.
Your main challenge in using this method will be to select your initial estimate for 'fzero' in such a way that it does not at some point set x7 to a negative quantity. That would produce a complex number. If necessary you could use as the input variable for Mattsfunc a quantity s where x7 = s^4 and the square root would be s^2, so there is no possibility of getting a complex result for sqrt(x7).
5 Comments
MATT
on 30 Aug 2014
Roger Stafford
on 30 Aug 2014
1. To "hand" means to write fzero(@Mattsfunc,x70) where the function "Mattsfunc" is so placed that it can also receive the parameters, n, c, and w, and where x70 is your initial estimate of the value of x7 which will hopefully give you a solution.
2. The matlab function 'fzero' will attempt to adjust x7 in such a way that F = x5^2*x3*x4-w becomes zero. That is what its job is. When (and if) it succeeds, you will have found a solution to all seven equations.
3. As I defined the function Mattsfunc, it is indeed a function of just the one variable x7. For each possible value of x7 there will be a unique value for F. x7 determines x5 and x6, they in turn determine x4, then x1, x2, x3, and finally F.
I would advise you to read the documentation for 'fzero' very carefully so that you understand how it works.
One observation: You will find that your system typically has either two different solutions for x7 or none for a given set of parameters, n, c, and w. In the first case you will need to experiment with differing initial estimates to enable fzero to find both solutions.
Roger Stafford
on 30 Aug 2014
Edited: Roger Stafford
on 30 Aug 2014
Just to reassure you, Matt, I obtained the following results using 'fzero' with my ancient version of matlab.
n = 345;
c = 12.3;
w = 0.1;
x7 = fzero(@Mattsfunc,3);
and got
x7 = 4.138413298429636
When I called
x7 = fzero(@Mattsfunc,7);
I received
x7 = 8.506854852971340
which demonstrates that your equations can have two solutions for the same set of parameters. In each case when the corresponding x1, x2, ..., x6 values are computed in accordance with the equations I gave in the function, both sets proved to be very accurate solutions to your seven equations.
I did have trouble with negative values for x7 where 'fzero' used an initial estimate of only +2. For that reason I would recommend using the interval form of the estimate so as to prevent negative values.
It helps greatly to make a plot of F = Mattsfunc(x7) as a function of x7. The curve is very smooth. It starts at F = -w for x7 = 0, rises to a peak, and finally converges asymptotically back down toward -w as x7 increases. It hopefully crosses zero at a couple of points in between. Such a curve can make it very much easier to make more accurate initial estimates. It can also indicate when there will be no solutions possible.
Roger Stafford
on 10 Sep 2014
There is a way of reducing seven (all but 1d.) of this new set of eight equations to a dependency on a single unknown, namely x8. However, when I generated several random sets of positive n, c, w, and co parameters and applied this method, no solutions could be found for any of them. The function that was derived from x8, namely x1+2*x2+x5-x3-2*x4, never crossed down past a zero value for positive values of x8, and that final equation 1d. could therefore not be satisfied.
I suspect that these new equations you have devised have no positive solutions for any set of positive parameters. However, if you are still interested in the method I can explain it to you. Perhaps you can suggest some likely ranges for the four parameters. Those I generated were all in the range from 0 to 1.
Roger Stafford
on 11 Sep 2014
I was overly pessimistic about finding solutions to your most recent set of equations. I had been testing the wrong collection of n, c, w, and co parameters. With appropriate values for them there are indeed solutions. Let's relabel your equations so the labels are all different:
1a. x1/x5 = n
1b. x6*x2/(x7*x5)^2 = c
1c. x7^2*x3*x5 = w
1d. x6*x4*(x7*x5)^2 = co
1e. x1+2*x2+x5-x3-2*x4 = 0
1f. x8 = 1/2*(x1+4*x2+x5-x3-4*x4)
1g. x7 = 10^(-x8^(1/2)/(1+x8^(1/2))+0.2*x8)
1h. x6 = 10^(-4*x8^(1/2)/(1+x8^(1/2))+0.8*x8)
The first step is to realize that, given equation 1e., then equation 1f. can be rewritten as:
1F. x8 = x2-x4
so we hereby replace 1f. by 1F. The solution space will remain the same.
If we choose some random value for x8, then x6 and x7 can be computed from 1h. and 1g. By multiplying 1b. and 1d. equations we get
x6^2*x2*x4 = c*co
and using 1F. gives us
(x4+x8)*x4 = c*co/y6^2
This is a quadratic equation in x4 and, assuming c and co are positive, there are two real roots, one positive and the other negative. We reject the negative one and get
x4 = (-x8+sqrt(x^2+c*co/y6^2)/2
x2 = x8 + x4 = (x8+sqrt(x^2+c*co/y6^2)/2
For more accurate computation the expression for x4 can be rewritten as the equivalent
x4 = 2*c*co/x6^2/(x8+sqrt(x8^2+4*c*co/x6^2))
by multiplying numerator and denominator by x8+sqrt(x8^2+4*c*co/x6^2).
Finally, using equations 1b., 1c., and 1a. we can find x5, x3, and x1.
Thus, starting with an arbitrary x8, we have derived all other x values and all equations will be satisfied except
1e. x1+2*x2+x5-x3-2*x4 = 0
Calculating x1+2*x2+x5-x3-2*x4 from x8 is the task that a function handed to 'fzero' can perform to search for the proper value of x8 that can make this expression zero and satisfy 1e. The following code carries out the above operations starting with given values of x8.
x7 = 10.^(-x8.^(1/2)./(1+x8.^(1/2))+0.2*x8);
x6 = 10.^(-4*x8.^(1/2)./(1+x8.^(1/2))+0.8*x8);
x4 = 2*c*co./x6.^2./(x8+sqrt(x8.^2+4*c*co./x6.^2));
x2 = (x8+sqrt(x8.^2+4*c*co./x6.^2))/2;
x5 = sqrt(co./(x6.*x4.*x7.^2));
x1 = n*x5;
x3 = w./(x7.^2.*x5);
y = x1+2*x2+x5-x3-2*x4;
The value of x8 is to be adjusted by 'fzero' until the value of y is zero.
I succeeded in obtaining a solution using these parameter values
n = 2;
c = 5;
w = 100;
co = 80;
and there should be quite an extensive set of possible parameter values that will also have solutions. However, as I found out earlier, not all values will succeed.
5 Comments
Roger Stafford
on 11 Sep 2014
That change does certainly invalidate the most recent method I presented. Your problem can still be solved by starting in with the single variable x8, and if it is done just right, you can reduce things to the solution of a certain cubic equation in x5, (not a quartic,) for which you can make use of matlab's 'roots' function. The difficulty there is that for certain combinations of parameters and trial values of x8, you may get more than one real positive root for x5, which makes coding more complicated.
At this point you might want to consider making use of matlab's 'fsolve' function using x8 and x5 as your two unknowns to be solved for. As you can easily see, your eight equations can readily be reduced to these two unknowns with two equations to solve. However, making initial estimates with this method is more difficult because making plots showing approximately where the roots lie using graphical techniques is much more complicated.
Roger Stafford
on 12 Sep 2014
Matt, I take back what I said about the supposed difficulty with the cubic equation. There is no such difficulty. Assuming that all four parameters are positive and that x8 is chosen positive, then it can be shown that there will always be exactly one and only one positive root for the cubic equation, so 'roots' will be easy to use. Do you want me to show you the details of this method?
MATT
on 7 Oct 2014
2 Comments
Sean de Wolski
on 17 Oct 2014
You could do this with fminbnd but you would have to square the output z before returning it in order to force 0 to be a minimum.
Roger Stafford
on 17 Oct 2014
Rather than placing a constraint on 'myfun' I would suggest calling on 'fzero' with a two-element "estimate" - that is, an interval in the x8 range. See the documentation at
http://www.mathworks.com/help/matlab/ref/fzero.html
which says: "x0 ... Initial value, specified as a real scalar or a 2-element real vector ... 2-element vector — fzero checks that fun(x0(1)) and fun(x0(2)) have opposite signs, and errors if they do not. It then iteratively shrinks the interval where fun changes sign to reach a solution."
If you have difficulties choosing this interval, you can use 'myfun' to make a plot of its 'z' output versus its 'x8' input over a wide range to see roughly where (or whether) z crosses zero. (Output 'z' is the quantity in eq. 5. that you are attempting to bring to zero.)
I noticed that for one set of n, c, w, and co values, z continued increasing above zero as x8 increases from 0 until it reached some maximum, and then decreased down past zero at higher values of x8, so you would want to have your interval chosen in that descending part on either side of the crossing point.
MATT
on 5 Nov 2014
Categories
Find more on Special Values in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
