solving a non-linear problem

9 vues (au cours des 30 derniers jours)
MATT
MATT le 16 Août 2014
Réponse apportée : MATT le 5 Nov 2014
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.

Connectez-vous pour répondre à cette question.

Réponse acceptée

Roger Stafford
Roger Stafford le 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.)

Plus de réponses (6)

Roger Stafford
Roger Stafford le 16 Août 2014
Modifié(e) : Roger Stafford le 16 Août 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.
  1 commentaire
MATT
MATT le 16 Août 2014
Hello, I am sorry, you are perfectly right I corrected the problem above. I forgot one equation. Thanks for reporting it.

Connectez-vous pour commenter.

Roger Stafford
Roger Stafford le 17 Août 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 commentaires
Afficher 1 commentaire plus ancien
Roger Stafford
Roger Stafford le 19 Août 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.
MATT
MATT le 29 Août 2014
Modifié(e) : MATT le 29 Août 2014
Hello Roger,
Thank you for your help. I have tried to apply the same method with the below system without success so far. I would like to use fzero to retrieve x4 first, and then use it to retrieve all other unknowns. Is it possible (could you also show me the development as you did before?)?
If it is not possible to use the fzero function for the present system, could you indicate me another way to solve the system with mathlab?
1a. x1/x4=n
1b. x5*x2/(x6*x4)^2=c
1c. x5^2*x3*x4=w
1d. x1+2*x2+x4-x3=0
1e. x7=1/2*(x1+4*x2+x3+x4)
1f. x5=10^(-x7^(1/2)/(1+x7^(1/2))+0.2*x7)
1g. x6=10^(-4*x7^(1/2)/(1+x7^(1/2))+0.8*x7)
Thanks so much,
Matthieu

Connectez-vous pour commenter.

Roger Stafford
Roger Stafford le 29 Août 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 commentaires
Afficher 3 commentaires plus anciens
MATT
MATT le 10 Sep 2014
Modifié(e) : MATT le 11 Sep 2014
Hello Roger, Thanks for all of this. It works perfectly and I have been able to solve many such systems. Now, I want to add one more layer of complexity with the following system
1a. x1/x5=n
1b. x6*x2/(x7*x5)^2=c
1c. x7^2*x3*x5=w
1c. x6*x4*(x7*x5)^2=co
1d. x1+2*x2+x5-x3-2*x4=0
1e. x8=1/2*(x1+4*x2+x5+x3+4*x4)
1f. x7=10^(-x8^(1/2)/(1+x8^(1/2))+0.2*x8)
1g. x6=10^(-4*x8^(1/2)/(1+x8^(1/2))+0.8*x8)
I have tried to proceed the same way you did but in this case it seems I fall on a 4th order polynomialto solve for x5, and I guess I need a guess for it as well. So from my tries, I suspect this system may need two guesses, one for x8 initially, and then for x5. Once again, all solutions and all parameters must be real and positive. If you still have energy, could you show me how to proceed with the steps. I then use these to redo the hand-calculation myself by adding many more equations in the system, all of similar forms than the one above that are representative. Best, Matthieu
Roger Stafford
Roger Stafford le 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.

Connectez-vous pour commenter.

Roger Stafford
Roger Stafford le 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 commentaires
Afficher 3 commentaires plus anciens
Roger Stafford
Roger Stafford le 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
MATT le 7 Oct 2014
Hello Roger, Sure, it would be great! I'd love to see the method for this case as well so that I can then reproduce it with other systems of similar forms. Also, I would need to fit a series of data points using an implicit non-linear function of the form f(x,y)=0. How can I do that? Thanks a lot Roger, Matthieu

Connectez-vous pour commenter.

MATT
MATT le 17 Oct 2014
Modifié(e) : MATT le 17 Oct 2014
|Hi Roger,
Beautiful. Is there a way to modify the "myfun" script in order to, rather than making an initial guess, stipulate a range in which the solution lies (which I can easily determine pretty accuratly) for x8 (between 0 and 3). Can we bracket zones of sign change and then converge to the solution?
Thanks a lot,
Matthieu|
  2 commentaires
Sean de Wolski
Sean de Wolski le 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
Roger Stafford le 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.

Connectez-vous pour commenter.

MATT
MATT le 5 Nov 2014
Dear Roger, I have the data points below, defined by the variables X1,X2,X3,X4:
X1 X2 X3 X4
3.5E-4 0 0.99965 1
4.08935E-4 0.00138 0.99821 0.85588
0.00114 0.03473 0.96413 0.30777
6.43932E-4 0.00989 0.98947 0.54354
0.00145 0.04977 0.94878 0.24193
0.0019 0.08999 0.90811 0.18414
0.00268 0.14986 0.84745 0.13042
0.00268 0.14998 0.84734 0.13056
0.00465 0.2488 0.74655 0.07521
I need to fit these data points using a function of the type
a*b*ln(X4)=W1*X2*(1-X1)+W2*X3*(1-X1)-W3*X2*X3+W4*X2*X3*(1-2*X1)
where a,b are constants, and W1,W2,W3,W4 are adjustable parameters. Could you indicate me how to do that? Best, Matthieu

Catégories

En savoir plus sur Interpolation dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by