solving a non-linear problem
3 views (last 30 days)
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.
0 Comments
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.)
0 Comments
More Answers (6)
Roger Stafford
on 16 Aug 2014
Edited: Roger Stafford
on 16 Aug 2014
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
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 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 17 Oct 2014
Edited: MATT
on 17 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.
See Also
Categories
Find more on Interpolation 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!