Thread Subject: Fsolve

Hi all,

i am trying to solve a system of linear/nonlinear equations using fsolve command in matlab. But i am getting some undesirable results.

First of all, kindly look at the code.

function F=myfun(X,c1)

a=X(1);
b=X(2);
c=X(3);
d=X(4);
e=X(5);
f=X(6);
g=X(7);
p=X(8);

F(1)=a+b+d+g-c1(4);
F(2)=c1(6)+2*c1(7)+c1(8)-2*a-b-f;
F(3)=c1(5)+2*c1(8)-2*c-4*d-2*f;
F(4)=3.76*c1(7)-e;
F(5)=c1(1)*a*p-(b*b);
F(6)=c1(2)*f*p-(b*c);
F(7)=c1(3)*c*c-(d*p);
F(8)=a+b+c+d+e+f-p;

as u can see from above function script, i have 8 unknowns (i.e a,b,c,d,e,f,g & p).
c1(1) to c1(8) are constant (having value greater then 0)

the fsolve command which i am using to solve the above function is,

fsolve(@(X) myfun_1(X,[q r s m x y z k]),[1 1 1 1 1 1 1 1],options);

where,
myfun_1 is function file name
q r s m x y z k are constant c1(1) to c1(8) (all greater than zero)
options = optimset('MaxFunEvals',3000,'MaxIter',20000);
[11111111] are initial guess

the problems i am facing with this code are,

1) most of the time the solution is not converged
         it gives me msgs like, change the initial guess ( which i changed thousands of time.)
         it ask me to increase the MaxFunEvals and iteration (it tried it to up to millions)

2) I am expecting values from this code between 0-5 for the variable (surely greater than zero) but it is giving me negative results.

so is some thing wrong with my code ?
is there any cure for it ?

kindly let me know any possible solution for this problem, if any other information is required then kindly let me know.

thanks and regard

have a nice day

bye.

rabbahs wrote:
> Hi all,
>
> i am trying to solve a system of linear/nonlinear equations using fsolve
> command in matlab. But i am getting some undesirable results.
>
> First of all, kindly look at the code.
>
> function F=myfun(X,c1)
>
> a=X(1);
> b=X(2);
> c=X(3);
> d=X(4);
> e=X(5);
> f=X(6);
> g=X(7);
> p=X(8);
>
> F(1)=a+b+d+g-c1(4);
> F(2)=c1(6)+2*c1(7)+c1(8)-2*a-b-f;
> F(3)=c1(5)+2*c1(8)-2*c-4*d-2*f;
> F(4)=3.76*c1(7)-e;
> F(5)=c1(1)*a*p-(b*b);
> F(6)=c1(2)*f*p-(b*c);
> F(7)=c1(3)*c*c-(d*p);
> F(8)=a+b+c+d+e+f-p;
>
> as u can see from above function script, i have 8 unknowns (i.e
> a,b,c,d,e,f,g & p).
> c1(1) to c1(8) are constant (having value greater then 0)

Exact solution:

Let P = RootOf(
10000*(4*c13 + 1)*(c11*c12^2*c13 + 4*c12^2*c13 - c11*c12 + c12^2 -
c11)*_Z^4 +

800*(25*c11*c12^2*c13*c16 + 169*c11*c12^2*c13*c17 - 188*c11*c12*c13*c17
+ 100*c12^2*c13*c16 + 1052*c12^2*c13*c17 - 25*c11*c12*c16 -
216*c11*c12*c17 + 25*c12^2*c16 + 263*c12^2*c17 - 25*c11*c16 -
169*c11*c17)*_Z^3 -

4*(625*c11*c12^2*c13*c15^2 + 2500*c11*c12^2*c13*c15*c16 +
16900*c11*c12^2*c13*c15*c17 + 2500*c11*c12^2*c13*c15*c18 +
5000*c11*c12^2*c13*c16*c18 + 33800*c11*c12^2*c13*c17*c18 +
2500*c11*c12^2*c13*c18^2 + 5000*c12^2*c13*c15^2 +
10000*c12^2*c13*c15*c16 + 105200*c12^2*c13*c15*c17 +
20000*c12^2*c13*c15*c18 + 20000*c12^2*c13*c16*c18 +
210400*c12^2*c13*c17*c18 + 20000*c12^2*c13*c18^2 - 625*c11*c12*c15^2 -
1250*c11*c12*c15*c16 - 13150*c11*c12*c15*c17 - 2500*c11*c12*c15*c18 +
2500*c11*c12*c16^2 + 43200*c11*c12*c16*c17 - 2500*c11*c12*c16*c18 +
177788*c11*c12*c17^2 - 26300*c11*c12*c17*c18 - 2500*c11*c12*c18^2 +
1250*c12^2*c15^2 + 2500*c12^2*c15*c16 + 26300*c12^2*c15*c17 +
5000*c12^2*c15*c18 - 2500*c12^2*c16^2 - 52600*c12^2*c16*c17 +
5000*c12^2*c16*c18 - 276676*c12^2*c17^2 + 52600*c12^2*c17*c18 +
5000*c12^2*c18^2 - 625*c11*c15^2 - 9400*c11*c15*c17 - 2500*c11*c15*c18 +
2500*c11*c16^2 + 33800*c11*c16*c17 + 78900*c11*c17^2 - 18800*c11*c17*c18
- 2500*c11*c18^2)*_Z^2 +

4*c12*(c15+2*c18)*(25*c15 + 50*c16 + 526*c17 + 50*c18)*(25*c11*c16 +
169*c11*c17 - 50*c12*c16 - 526*c12*c17)*_Z +
625*c12^2*c15^4 + 2500*c12^2*c15^3*c16 + 26300*c12^2*c15^3*c17 +
5000*c12^2*c15^3*c18 + 2500*c12^2*c15^2*c16^2 +
52600*c12^2*c15^2*c16*c17 + 15000*c12^2*c15^2*c16*c18 +
276676*c12^2*c15^2*c17^2 + 157800*c12^2*c15^2*c17*c18 +
15000*c12^2*c15^2*c18^2 + 10000*c12^2*c15*c16^2*c18 +
210400*c12^2*c15*c16*c17*c18 + 30000*c12^2*c15*c16*c18^2 +
1106704*c12^2*c15*c17^2*c18 + 315600*c12^2*c15*c17*c18^2 +
20000*c12^2*c15*c18^3 + 10000*c12^2*c16^2*c18^2 +
210400*c12^2*c16*c17*c18^2 + 20000*c12^2*c16*c18^3 +
1106704*c12^2*c17^2*c18^2 + 210400*c12^2*c17*c18^3 + 10000*c12^2*c18^4 )

The RootOf() operator means to find the values of _Z that make the
expression 0. Notice that this is a quartic expression in _Z and thus
has exact solutions -- but printing out the exact solutions would be
rather long! Easier to substitute in your known constants and use
Matlab's roots() operator to get the roots -- some of which may well be
complex. And you might need the complex roots in order to proceed further.

Once P are determined,

Let Q = 50*(2*c12*c13 - 1)*P - 25*c15 - 50*c16 - 526*c17-50*c18
Let R = 2*(c12 + 2)*P - c12*c15 - 2*c12*c18

Then:

a = (2500*(4*c13+1)*(c12^2*c13 - c12 - 1)*P^3 + (5000*c12^2*c13*c16 +
33800*c12^2*c13*c17 - 37600*c12*c13*c17 - 5000*c12*c16 - 43200*c12*c17 -
5000*c16 - 33800*c17)*P^2 + (-625*c12^2*c13*c15^2 -
2500*c12^2*c13*c15*c16 - 16900*c12^2*c13*c15*c17 -
2500*c12^2*c13*c15*c18 - 5000*c12^2*c13*c16*c18 -
33800*c12^2*c13*c17*c18 - 2500*c12^2*c13*c18^2 + 625*c12*c15^2 +
1250*c12*c15*c16 + 13150*c12*c15*c17 + 2500*c12*c15*c18 - 2500*c12*c16^2
- 43200*c12*c16*c17 + 2500*c12*c16*c18 - 177788*c12*c17^2 +
26300*c12*c17*c18 + 2500*c12*c18^2 + 625*c15^2 + 9400*c15*c17 +
2500*c15*c18 - 2500*c16^2 - 33800*c16*c17 - 78900*c17^2 + 18800*c17*c18
+ 2500*c18^2)*P + 625*c12*c15^2*c16 + 4225*c12*c15^2*c17 +
1250*c12*c15*c16^2 + 21600*c12*c15*c16*c17 + 2500*c12*c15*c16*c18 +
88894*c12*c15*c17^2 + 16900*c12*c15*c17*c18 + 2500*c12*c16^2*c18 +
43200*c12*c16*c17*c18 + 2500*c12*c16*c18^2 + 177788*c12*c17^2*c18 +
16900*c12*c17*c18^2) / (25*Q*R)

b = ( -100*c12*(4*c13 + 1)*P^2 - 4*c12*(263*c17 + 25*c16)*P +
c12*(c15+2*c18)*(25*c15+50*c16 + 526*c17 + 50*c18) ) / (50*R)

c = P

d = -25*c13*P*R/Q

e = 94/25*c17

f = 1/2* (100*(4*c13 + 1)*P^2 + 4*(25*c16 + 263*c17)*P - 25*c15^2 -
50*c15*c16 - 526*c15*c17 - 100*c15*c18 - 100*c16*c18 - 1052*c17*c18 -
100*c18^2) / ((100*c12*c13 - 50)*P - 25*c15 - 50*c16 - 526*c17 - 50*c18)

g = ( 5000*(4*c12^2*c13^2 + 2*c12^2*c13 + 4*c12*c13 + 8*c13 + 1)*P^3 +
100*(100*c12^2*c13*c14 - 50*c12^2*c13*c15 + 376*c12^2*c13*c17 -
100*c12^2*c13*c18 + 200*c12*c13*c14 - 200*c12*c13*c15 - 200*c12*c13*c16
- 1352*c12*c13*c17 - 400*c12*c13*c18 - 50*c12*c14-25*c12*c15 -
188*c12*c17 - 50*c12*c18 - 100*c14+100*c16 + 676*c17)*P^2 -
2*(2500*c12^2*c13*c14*c15 + 5000*c12^2*c13*c14*c18 +
9400*c12^2*c13*c15*c17 + 18800*c12^2*c13*c17*c18 + 2500*c12*c14*c16 +
26300*c12*c14*c17 + 1250*c12*c15*c16 + 13150*c12*c15*c17 +
9400*c12*c16*c17 + 2500*c12*c16*c18 + 98888*c12*c17^2 +
26300*c12*c17*c18 + 2500*c14*c15 + 5000*c14*c16 + 52600*c14*c17 +
5000*c14*c18 + 625*c15^2 + 9400*c15*c17 + 2500*c15*c18 - 2500*c16^2 -
33800*c16*c17 - 78900*c17^2 + 18800*c17*c18 + 2500*c18^2)*P +
9400*c12*c15*c16*c17 + 71400*c12*c15*c17*c18 + 5000*c12*c15*c16*c18 +
18800*c12*c16*c17*c18 + 625*c12*c15^3 + 5000*c12*c18^3 +
98888*c12*c15*c17^2 + 17850*c12*c15^2*c17 + 3750*c12*c15^2*c18 +
7500*c12*c15*c18^2 + 1250*c12*c15^2*c16 + 5000*c12*c16*c18^2 +
197776*c12*c17^2*c18 + 71400*c12*c17*c18^2 + 26300*c12*c14*c15*c17 +
5000*c12*c14*c15*c18 + 5000*c12*c14*c16*c18 + 52600*c12*c14*c17*c18 +
1250*c12*c14*c15^2 + 5000*c12*c14*c18^2 + 2500*c12*c14*c15*c16 ) / (50*Q*R)

p = -1/25*P*Q / ( (2*c12 + 4)*P - c12*c15 - 2*c12*c18 )

If you examine b, d, and p, you can see that it wouldn't take much for
the values to come out negative.


Substituting in the random constants [5/2, 9/5, 33/10, 31/10, 7/5,
17/10, 17/5, 23/10], I get the results:

[a = 1.267103339, b = 9.012979781, c = 1.806899943, d = .4201432576, e =
12.78400000, f = .3528135419, g = -7.600226377, p = 25.64393986]

[a = 18.20425353, b = -21.83226313, c = 2.310915533, d = 1.682664193, e
= 12.78400000, f = -2.676243918, g = 5.045345414, p = 10.47332620]

[a = 14.20536180, b = -18.69989687, c = -2.074911341, d = 1.442869037, e
= 12.78400000, f = 2.189173266, g = 6.151666031, p = 9.846595894]

[a = 1.749657476, b = 9.856441847, c = -5.904572379, d = 5.180164589, e
= 12.78400000, f = -1.455756799, g = -13.68626391, p = 22.20993473]

Notice that no matter which root of the quartic was chosen, at least one
of the variables came out negative. I do not know at the moment if that
is *necessarily* the case -- my computation engine has pretty much
locked itself up processing the testing. I can say this: g is a function
of c1(8) to the 12th power, involving quantities ranging from 10^(-10)
to 10^54. With that range of values, order of operations would become
crucial as round-off error is likely to *very* high. Scanning through
the expression, if you were to calculate it naively in double precision,
I can see that you would be lucky to get a single bit of the solution
right. This high reliance upon the exact value of c18 is hidden in the
way I rewrote the solution, above: P, the RootOf, involves up to c18^4,
and g depends upon P^3 and thus upon (c18^4)^3 = c18^12 .

Mr. Walter,

I really appreciate your effort to put your time to solve my problem. thanks

what i get from you reply is that you find out the exact solution for the variable using Rootof( command in Matlab (i hope this is the case).

what are these random constant ? c11 c12 ...


as you find out that one of the value among variable always comes negative (does it mean that there is some problem in the equations it self, thats y it is giving negative values ?)

how can i write the nonlinear equations in AX=B form.

kindly refer me to some tutorial for fsolve (other then Matlab tutorial)

again thanks

have a nice day
bye

rabbahs wrote:

 > how can i write the nonlinear equations in AX=B form.

You can't. Some of your variables appear both squared and to the first
power, so you cannot handle the situation by variable substitution.


> what i get from you reply is that you find out the exact solution for
> the variable using Rootof( command in Matlab (i hope this is the case).

The Matlab equivalent of Maple's RootOf is roots()


> what are these random constant ? c11 c12 ...

c11 is what you wrote as c1(1)
c12 is what you wrote as c1(2)
and so on. I needed to rename during my processing because otherwise my
symbolic engine would have thought that I was indicating function calls.


> as you find out that one of the value among variable always comes
> negative (does it mean that there is some problem in the equations it
> self, thats y it is giving negative values ?)

I discovered a moment ago that I had a typo when I input your equations
originally, so the solution I gave is not valid!!

I'm working towards a more compact solution now... or at least a more
correct solution.

Mr. Walter,

I really appreciate your effort to put your time to solve my problem. thanks

what i get from you reply is that you find out the exact solution for the variable using Rootof( command in Matlab (i hope this is the case).

what are these random constant ? c11 c12 ...


as you find out that one of the value among variable always comes negative (does it mean that there is some problem in the equations it self, thats y it is giving negative values ?)

how can i write the nonlinear equations in AX=B form.

kindly refer me to some tutorial for fsolve (other then Matlab tutorial)

again thanks

have a nice day
bye

Thanks Mr. Walter for your time.

I m looking fwd to see the solution.

If you want then i can send you my script file ? (it might help you to figure out the problem easily).

have a nice day

bye

rabbahs wrote:
> Thanks Mr. Walter for your time.
> I m looking fwd to see the solution.
>
> If you want then i can send you my script file ? (it might help you to
> figure out the problem easily).

I have determined that it is not _necessarily_ the case that some
variables must come out negative. It depends on the initial constants.

For the constants
c1(1) = 5/2
c1(2) = 9/5
c1(3) = 33/10
c1(4) = 10
c1(5) = 7/5
c1(6) = 17/10
c1(7) = 17/5
c1(8) = 23/10

I reach the solution

a = 1.103817340
b = 8.258157698
d = .4327670834
c = 1.800258210
e = 12.784
f = .33420762
g = .205257879
p = 24.71320795

The solution process involved reducing it down to a problem in 3
expressions, with the other 5 reducing to 0 by linear combinations of
the 3 remaining variables. Those expressions (which have to be solved to
equal 0) are:

c11*a*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b^2,

c12*(c16 + 2*c17 + c18 - 2*a - b)*(a + b + 1/2*c15 - d + 94/25*c17 +
c18) - b*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a+b),

c13*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b)^2 - d*(a + b + 1/2*c15 - d +
94/25*c17 + c18)

As before, c18 is what you wrote as the constant c1(8) and so on.

The solution to this trio of equations involves finding the roots of a
quartic. The generalized solution is pretty long and messy, so it took
me numerous false starts to find a tractable approach. What I did was
generate random values for all of the constants *except* c1(4), and
substituted those values in; c1(4) is not present in the trio of
equations, so it is possible to find values for the four tuples of
solutions for a, b, and d. Four tuples because a quartic has four
solutions. With the particular random values I had, two of the tuples
had negative values for b, so those two tuples could be eliminated
immediately. The other two tuples were all positive for a, b, and d. I
then back-substituted a, b, and d in to find the values for c, e, f, and
p, and to find the expression for g, which is c14 - a - b - d . One of
the two tuples gave negative values for c and f; the other gave all
positive values for the variables, with g = c14 minus a constant. From
there it was just a matter of selecting any c14 larger than the constant
to arrive at a set of all-positive constants that will produce a set of
all-positive variables.

This is proof by construction: I constructed a set of constants that was
solvable for solutions in the required range. As such, it is only a
proof that desirable solutions exist for _some_ combinations of
constants: other combinations of constants will have no valid solutions.


The values that can be deduced linearly from a, b, and d are:

c = 1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b
e = 94/25*c17
f = c16 + 2*c17 + c18 - 2*a - b
g = c14 - a - b - d
p = a + b + c + d + e + f

p can also be expressed as p = a + b - d + c15/2 + 94/25*c17 + c18
by combining the other equations.

By examination, we can see that e is the only variable which is
guaranteed to come out positive based upon the constants; everything
else depends upon a combination of the values of the constants and the
solutions found to the 3 non-linear equations. And make sure to check
all four roots of the quartic, as they might not meet your constraints.

When you are calculating the roots of the quartic, it is common that
most or all of the roots you get out will be imaginary numbers. If the
imaginary components are on the order of 10^(-6) or less, then those are
being introduced by round-off error and you should just use real() to
extract the real portions. But you do have to be careful, because with
some values of the constants, you can get roots whose imaginary portions
are on the order of .05 to -21 (values observed in practice), and those
values are much too large for round off: those are true imaginary
numbers and taking the real portion of them will seriously mislead you.

Walter Roberson <roberson@hushmail.com> wrote in message <ho4itl$msf$1@canopus.cc.umanitoba.ca>...
> rabbahs wrote:
> > Thanks Mr. Walter for your time.
> > I m looking fwd to see the solution.
> >
> > If you want then i can send you my script file ? (it might help you to
> > figure out the problem easily).
>
> I have determined that it is not _necessarily_ the case that some
> variables must come out negative. It depends on the initial constants.
>
> For the constants
> c1(1) = 5/2
> c1(2) = 9/5
> c1(3) = 33/10
> c1(4) = 10
> c1(5) = 7/5
> c1(6) = 17/10
> c1(7) = 17/5
> c1(8) = 23/10
>
> I reach the solution
>
> a = 1.103817340
> b = 8.258157698
> d = .4327670834
> c = 1.800258210
> e = 12.784
> f = .33420762
> g = .205257879
> p = 24.71320795
>
> The solution process involved reducing it down to a problem in 3
> expressions, with the other 5 reducing to 0 by linear combinations of
> the 3 remaining variables. Those expressions (which have to be solved to
> equal 0) are:
>
> c11*a*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b^2,
>
> c12*(c16 + 2*c17 + c18 - 2*a - b)*(a + b + 1/2*c15 - d + 94/25*c17 +
> c18) - b*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a+b),
>
> c13*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b)^2 - d*(a + b + 1/2*c15 - d +
> 94/25*c17 + c18)
>
> As before, c18 is what you wrote as the constant c1(8) and so on.
>
> The solution to this trio of equations involves finding the roots of a
> quartic. The generalized solution is pretty long and messy, so it took
> me numerous false starts to find a tractable approach. What I did was
> generate random values for all of the constants *except* c1(4), and
> substituted those values in; c1(4) is not present in the trio of
> equations, so it is possible to find values for the four tuples of
> solutions for a, b, and d. Four tuples because a quartic has four
> solutions. With the particular random values I had, two of the tuples
> had negative values for b, so those two tuples could be eliminated
> immediately. The other two tuples were all positive for a, b, and d. I
> then back-substituted a, b, and d in to find the values for c, e, f, and
> p, and to find the expression for g, which is c14 - a - b - d . One of
> the two tuples gave negative values for c and f; the other gave all
> positive values for the variables, with g = c14 minus a constant. From
> there it was just a matter of selecting any c14 larger than the constant
> to arrive at a set of all-positive constants that will produce a set of
> all-positive variables.
>
> This is proof by construction: I constructed a set of constants that was
> solvable for solutions in the required range. As such, it is only a
> proof that desirable solutions exist for _some_ combinations of
> constants: other combinations of constants will have no valid solutions.
>
>

Thanks Mr. Walter,

I m working on the problem, will let you know the results as soon as i done with it.

thanks for all the help !!!

If i have any query, i will let you know.
> The values that can be deduced linearly from a, b, and d are:
>
> c = 1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b
> e = 94/25*c17
> f = c16 + 2*c17 + c18 - 2*a - b
> g = c14 - a - b - d
> p = a + b + c + d + e + f
>
> p can also be expressed as p = a + b - d + c15/2 + 94/25*c17 + c18
> by combining the other equations.
>
> By examination, we can see that e is the only variable which is
> guaranteed to come out positive based upon the constants; everything
> else depends upon a combination of the values of the constants and the
> solutions found to the 3 non-linear equations. And make sure to check
> all four roots of the quartic, as they might not meet your constraints.
>
> When you are calculating the roots of the quartic, it is common that
> most or all of the roots you get out will be imaginary numbers. If the
> imaginary components are on the order of 10^(-6) or less, then those are
> being introduced by round-off error and you should just use real() to
> extract the real portions. But you do have to be careful, because with
> some values of the constants, you can get roots whose imaginary portions
> are on the order of .05 to -21 (values observed in practice), and those
> values are much too large for round off: those are true imaginary
> numbers and taking the real portion of them will seriously mislead you.

rabbahs wrote:
> Walter Roberson <roberson@hushmail.com> wrote in message
> <ho4itl$msf$1@canopus.cc.umanitoba.ca>...

>> The values that can be deduced linearly from a, b, and d are:
>>
>> c = 1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b
>> e = 94/25*c17
>> f = c16 + 2*c17 + c18 - 2*a - b
>> g = c14 - a - b - d
>> p = a + b + c + d + e + f
>>
>> p can also be expressed as p = a + b - d + c15/2 + 94/25*c17 + c18
>> by combining the other equations.

If you examine the expressions for c, g, and p, you can see that d
always appears in the negative for them. Thus those expressions have the
greatest chance of being positive if d is as small as possible. To
explore that possibility, evaluate the system of three non-linear
equations at d = 0, leaving a system of three equations in two
variables. The second of the equations is the most complex of them, so
solve the first against the third. The result is a relatively simple
quadratic which lead to "reasonable" expressions for a and b -- indeed,
one of the two solutions to the quadratic has b = 0. Choosing the other
root and substituting back in to the system of three equations, sure
enough the first and third go to 0, and the second goes to something
that is not all that long -- and which factors moderately well.

If you then check which constants that expression involve, and solve the
expression for each of them in turn (one at a time), then you find that
each in turn is the sum of negatives of some of the constants. Under the
assumption that the constants are positive, the implication is that the
constant being solved for would end up negative; and since this happens
for all the constants involved, there is no possibility to solve the
system to end up with all positive constants. From this, we deduce that
in fact d must be non-zero and large enough to interact with the other
constants to produce magnitudes comparable to the other constants.

I proceeded to experiment with the same solving method for d=1, but the
equations it spits out are horribly long, and have my symbolic
computation engine pretty much locked up trying to simplify them... the
fourth common subexpression is somewhere over 6000 lines long, and the
main expression and earlier common sub-expressions have scrolled off my
10000 line scrollback...

Walter Roberson <roberson@hushmail.com> wrote in message <ho6nln$bnj$1@canopus.cc.umanitoba.ca>...
> rabbahs wrote:
> > Walter Roberson <roberson@hushmail.com> wrote in message
> > <ho4itl$msf$1@canopus.cc.umanitoba.ca>...
>
> >> The values that can be deduced linearly from a, b, and d are:
> >>
> >> c = 1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b
> >> e = 94/25*c17
> >> f = c16 + 2*c17 + c18 - 2*a - b
> >> g = c14 - a - b - d
> >> p = a + b + c + d + e + f
> >>
> >> p can also be expressed as p = a + b - d + c15/2 + 94/25*c17 + c18
> >> by combining the other equations.
>
> If you examine the expressions for c, g, and p, you can see that d
> always appears in the negative for them. Thus those expressions have the
> greatest chance of being positive if d is as small as possible. To
> explore that possibility, evaluate the system of three non-linear
> equations at d = 0, leaving a system of three equations in two
> variables. The second of the equations is the most complex of them, so
> solve the first against the third. The result is a relatively simple
> quadratic which lead to "reasonable" expressions for a and b -- indeed,
> one of the two solutions to the quadratic has b = 0. Choosing the other
> root and substituting back in to the system of three equations, sure
> enough the first and third go to 0, and the second goes to something
> that is not all that long -- and which factors moderately well.
>
> If you then check which constants that expression involve, and solve the
> expression for each of them in turn (one at a time), then you find that
> each in turn is the sum of negatives of some of the constants. Under the
> assumption that the constants are positive, the implication is that the
> constant being solved for would end up negative; and since this happens
> for all the constants involved, there is no possibility to solve the
> system to end up with all positive constants. From this, we deduce that
> in fact d must be non-zero and large enough to interact with the other
> constants to produce magnitudes comparable to the other constants.
>
> I proceeded to experiment with the same solving method for d=1, but the
> equations it spits out are horribly long, and have my symbolic
> computation engine pretty much locked up trying to simplify them... the
> fourth common subexpression is somewhere over 6000 lines long, and the
> main expression and earlier common sub-expressions have scrolled off my
> 10000 line scrollback...



Hi Mr. Walter,

can u give me the expression to calculate the positive value of a,b and d ?

rabbahs wrote:

> can u give me the expression to calculate the positive value of a,b and d ?

Well, as I wrote before:

====
The solution process involved reducing it down to a problem in 3
expressions, with the other 5 reducing to 0 by linear combinations of
the 3 remaining variables. Those expressions (which have to be solved to
equal 0) are:

c11*a*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b^2,

c12*(c16 + 2*c17 + c18 - 2*a - b)*(a + b + 1/2*c15 - d + 94/25*c17 +
c18) - b*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a+b),

c13*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b)^2 - d*(a + b + 1/2*c15 - d +
94/25*c17 + c18)

As before, c18 is what you wrote as the constant c1(8) and so on.

The solution to this trio of equations involves finding the roots of a
quartic.
=====

So substitute your known constants into the three equations above, and
then proceed with traditional elimination. For example, d appears in the
first one multiplied only by a and a constant, so move the -b^2 to the
right hand side, divide by c11*a, and you then have a linear expression
on the left hand side that can be used to isolate d. That d value can
then be substituted into the two remaining equations. The second
expression will then become cubic in b but the third expression will
become quartic in b, so it would probably be easier to solve the second
expression first for a (make sure you consider all three roots), and
then use that to solve the third expression.


After substitution of the constants, the three expressions will be of
the general form

E1*a*(a+b+E2-d)/E3-b^2,

E4*(E5-2*a-b)*(a+b+E2-d)/E6-b*(E7-2*d+2*a+b),

E8*(E9-2*d+2*a+b)/E10-d*(a+b+E2-d)

where the E's are a combination of constants determined by the c1(:)
values. There is a general symbolic solution to this, but the general
solution has a repeated subexpression a quartic polynomial. In Matlab
notation, this quartic can be solved via

P = roots([E10*(E1+4*E3) * (E1*E6*E4+E1*E6^2-E3*E4^2),

   E1*(8*E8*E4*E3*E6 - 2*E10*E2*E4^2*E3 - 4*E10*E3*E6*E4*E5 -
     E1*E10*E6*E7*E4 - 4*E10*E6*E7*E4*E3 + E1*E8*E4^2 + 3*E1*E10*E2*E6*E4
     + 4*E10*E2*E4*E3*E6 + 2*E1*E8*E4*E6 - E1*E10*E6*E4*E5 +
     2*E1*E10*E2*E6^2 + 4*E8*E4^2*E3),

   E1*(E1*E10*E6*E7*E4*E5 + 4*E1*E8*E2*E4*E6 - 2*E1*E8*E4*E5*E6 +
     E10*E3*E4^2*E5^2 + 2*E1*E10*E2*E6^2*E7 + 4*E1*E8*E6^2*E7 -
     E1*E10*E6^2*E7^2 + 2*E1*E8*E6*E7*E4 - 4*E8*E9*E1*E6^2 -
     4*E8*E4^2*E3*E5 + 2*E1*E10*E2^2*E6*E4 - E8*E9*E1*E4^2 +
     2*E1*E8*E2*E4^2 + 2*E10*E2*E4^2*E3*E5 - 3*E1*E10*E2*E6*E4*E5 -
     2*E1*E8*E4^2*E5 - E1*E10*E2*E6*E7*E4 - 4*E8*E9*E1*E4*E6),

   E1^2*E4*E5*(-E8*E4*E5 - 2*E4*E8*E9 + 4*E8*E2*E4 + 2*E8*E6*E7 -
     4*E8*E9*E6 + 2*E10*E2^2*E6 - E10*E2*E6*E7 + 4*E8*E2*E6),

   -E8*E9*E1^2*E4^2*E5^2 + 2*E8*E1^2*E2*E4^2*E5^2 )

where the various E's are copied by inspection of the coefficients of
the three equations above after substitution of the constants.

Once you have the roots of this P value, substitute them each in turn into:

a = -P*E3*((-E4 - 2*E6)*P + E4*E5) /
  ((E1*E6 - 2*E3*E4)*P - E6*E7*E1 + 2*E6*E1*E2)

b = P

d = (-(E1 + 4*E3)*(E1*E6*E4 + E1*E6^2 - E3*E4^2)*P^3 - E1*(3*E1*E4*E6*E2
   - E1*E4*E5*E6 - 2*E2*E4^2*E3 - E1*E4*E6*E7 - 4*E3*E6*E4*E5 +
   4*E2*E6*E4*E3 + 2*E1*E6^2*E2 - 4*E6*E7*E4*E3)*P^2 - E1*(E3*E4^2*E5^2 -
   E1*E6^2*E7^2 - 3*E1*E4*E5*E6*E2 + 2*E1*E2*E6^2*E7 + 2*E1*E2^2*E4*E6 -
   E1*E2*E4*E6*E7 + E1*E4*E5*E6*E7 + 2*E2*E4^2*E5*E3)*P -
   E1^2*E2*E4*E5*E6*E7 + 2*E1^2*E2^2*E4*E5*E6) /
  (((-E4 - 2*E6)*P + E4*E5)*E1*((E1*E6 - 2*E3*E4)*P - E6*E7*E1 +
     2*E6*E1*E2)

and see if you can find a root set of P that makes a, b, and d all
non-negative. If none of the root-sets work, then there is no hope with
that combination of constants. (This set of expressions was found via
the Maple symbolic combination engine -- I don't try to do these things
by hand!)

Once you have an all-positive (a, b, d) tuple, substitute the respective
constants and a, b, d values into

c = 1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b
e = 94/25*c17
f = c16 + 2*c17 + c18 - 2*a - b
g = c14 - a - b - d
p = a + b - d + c15/2 + 94/25*c17 + c18

and see if you get all-positive results. If so, you are done. If not,
then go back and try the next root-set of the output of the roots()
call. If none of the four work, then you cannot solve the problem with
that combination of constants.

I do not recommend trying to solve this system via a minimizer: the
constants either have the right relationship for the exact solution, or
they do not have the right relationship and no amount of minimization
would fix that. The only time I might suggest a minimizer on this is if
the "constants" are actually intervals and you are hoping that somewhere
in the intervals is a solution such that all the variables are positive.


If you explore this system of equations further, you will find that in
general it has singularities. The singularities _tend_ to be in the
negative range, but with some combinations of values, the singularities
can occur in the positive ranges. For example, with the set of random
c1(:) constants I indicated before that I posted with, except with c18
left undetermined, I found that the graph was quite wobbly between 8.1
and 8.4 with multiple roots _and_ a singularity within a range of about
0.102. The singularity itself, Matlab cannot isolate, because a change
in the last bit causes a jump from about -21 to about +7: there might
actually be a zero in there, if you used used much higher precision...
but be off by .1 and the values jump between +/- 10^64 !

Walter Roberson <roberson@hushmail.com> wrote in message <ho7ejp$hro$1@canopus.cc.umanitoba.ca>...
> rabbahs wrote:
>
> > can u give me the expression to calculate the positive value of a,b and d ?
>
> Well, as I wrote before:
>
> ====
> The solution process involved reducing it down to a problem in 3
> expressions, with the other 5 reducing to 0 by linear combinations of
> the 3 remaining variables. Those expressions (which have to be solved to
> equal 0) are:
>
> c11*a*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b^2,
>
> c12*(c16 + 2*c17 + c18 - 2*a - b)*(a + b + 1/2*c15 - d + 94/25*c17 +
> c18) - b*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a+b),
>
> c13*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b)^2 - d*(a + b + 1/2*c15 - d +
> 94/25*c17 + c18)
>
> As before, c18 is what you wrote as the constant c1(8) and so on.
>
> The solution to this trio of equations involves finding the roots of a
> quartic.
> =====
>
> So substitute your known constants into the three equations above, and
> then proceed with traditional elimination. For example, d appears in the
> first one multiplied only by a and a constant, so move the -b^2 to the
> right hand side, divide by c11*a, and you then have a linear expression
> on the left hand side that can be used to isolate d. That d value can
> then be substituted into the two remaining equations. The second
> expression will then become cubic in b but the third expression will
> become quartic in b, so it would probably be easier to solve the second
> expression first for a (make sure you consider all three roots), and
> then use that to solve the third expression.
>
>
> After substitution of the constants, the three expressions will be of
> the general form
>
> E1*a*(a+b+E2-d)/E3-b^2,
>
> E4*(E5-2*a-b)*(a+b+E2-d)/E6-b*(E7-2*d+2*a+b),
>
> E8*(E9-2*d+2*a+b)/E10-d*(a+b+E2-d)
>
> where the E's are a combination of constants determined by the c1(:)
> values. There is a general symbolic solution to this, but the general
> solution has a repeated subexpression a quartic polynomial. In Matlab
> notation, this quartic can be solved via
>
> P = roots([E10*(E1+4*E3) * (E1*E6*E4+E1*E6^2-E3*E4^2),
>
> E1*(8*E8*E4*E3*E6 - 2*E10*E2*E4^2*E3 - 4*E10*E3*E6*E4*E5 -
> E1*E10*E6*E7*E4 - 4*E10*E6*E7*E4*E3 + E1*E8*E4^2 + 3*E1*E10*E2*E6*E4
> + 4*E10*E2*E4*E3*E6 + 2*E1*E8*E4*E6 - E1*E10*E6*E4*E5 +
> 2*E1*E10*E2*E6^2 + 4*E8*E4^2*E3),
>
> E1*(E1*E10*E6*E7*E4*E5 + 4*E1*E8*E2*E4*E6 - 2*E1*E8*E4*E5*E6 +
> E10*E3*E4^2*E5^2 + 2*E1*E10*E2*E6^2*E7 + 4*E1*E8*E6^2*E7 -
> E1*E10*E6^2*E7^2 + 2*E1*E8*E6*E7*E4 - 4*E8*E9*E1*E6^2 -
> 4*E8*E4^2*E3*E5 + 2*E1*E10*E2^2*E6*E4 - E8*E9*E1*E4^2 +
> 2*E1*E8*E2*E4^2 + 2*E10*E2*E4^2*E3*E5 - 3*E1*E10*E2*E6*E4*E5 -
> 2*E1*E8*E4^2*E5 - E1*E10*E2*E6*E7*E4 - 4*E8*E9*E1*E4*E6),
>
> E1^2*E4*E5*(-E8*E4*E5 - 2*E4*E8*E9 + 4*E8*E2*E4 + 2*E8*E6*E7 -
> 4*E8*E9*E6 + 2*E10*E2^2*E6 - E10*E2*E6*E7 + 4*E8*E2*E6),
>
> -E8*E9*E1^2*E4^2*E5^2 + 2*E8*E1^2*E2*E4^2*E5^2 )
>
> where the various E's are copied by inspection of the coefficients of
> the three equations above after substitution of the constants.
>
> Once you have the roots of this P value, substitute them each in turn into:
>
> a = -P*E3*((-E4 - 2*E6)*P + E4*E5) /
> ((E1*E6 - 2*E3*E4)*P - E6*E7*E1 + 2*E6*E1*E2)
>
> b = P
>
> d = (-(E1 + 4*E3)*(E1*E6*E4 + E1*E6^2 - E3*E4^2)*P^3 - E1*(3*E1*E4*E6*E2
> - E1*E4*E5*E6 - 2*E2*E4^2*E3 - E1*E4*E6*E7 - 4*E3*E6*E4*E5 +
> 4*E2*E6*E4*E3 + 2*E1*E6^2*E2 - 4*E6*E7*E4*E3)*P^2 - E1*(E3*E4^2*E5^2 -
> E1*E6^2*E7^2 - 3*E1*E4*E5*E6*E2 + 2*E1*E2*E6^2*E7 + 2*E1*E2^2*E4*E6 -
> E1*E2*E4*E6*E7 + E1*E4*E5*E6*E7 + 2*E2*E4^2*E5*E3)*P -
> E1^2*E2*E4*E5*E6*E7 + 2*E1^2*E2^2*E4*E5*E6) /
> (((-E4 - 2*E6)*P + E4*E5)*E1*((E1*E6 - 2*E3*E4)*P - E6*E7*E1 +
> 2*E6*E1*E2)
>
> and see if you can find a root set of P that makes a, b, and d all
> non-negative. If none of the root-sets work, then there is no hope with
> that combination of constants. (This set of expressions was found via
> the Maple symbolic combination engine -- I don't try to do these things
> by hand!)
>
> Once you have an all-positive (a, b, d) tuple, substitute the respective
> constants and a, b, d values into
>
> c = 1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b
> e = 94/25*c17
> f = c16 + 2*c17 + c18 - 2*a - b
> g = c14 - a - b - d
> p = a + b - d + c15/2 + 94/25*c17 + c18
>
> and see if you get all-positive results. If so, you are done. If not,
> then go back and try the next root-set of the output of the roots()
> call. If none of the four work, then you cannot solve the problem with
> that combination of constants.
>
> I do not recommend trying to solve this system via a minimizer: the
> constants either have the right relationship for the exact solution, or
> they do not have the right relationship and no amount of minimization
> would fix that. The only time I might suggest a minimizer on this is if
> the "constants" are actually intervals and you are hoping that somewhere
> in the intervals is a solution such that all the variables are positive.
>
>
> If you explore this system of equations further, you will find that in
> general it has singularities. The singularities _tend_ to be in the
> negative range, but with some combinations of values, the singularities
> can occur in the positive ranges. For example, with the set of random
> c1(:) constants I indicated before that I posted with, except with c18
> left undetermined, I found that the graph was quite wobbly between 8.1
> and 8.4 with multiple roots _and_ a singularity within a range of about
> 0.102. The singularity itself, Matlab cannot isolate, because a change
> in the last bit causes a jump from about -21 to about +7: there might
> actually be a zero in there, if you used used much higher precision...
> but be off by .1 and the values jump between +/- 10^64 !

Thanks Mr. Walter,

currently i m using fmincon command to solve this problem. I reduced the equations in to 3 set of equations (with variable a,b and 9, as u advise earlier) and that try to bound my solution with the constrains available with this function.

here are some details,

global c11 c12 c13 c14 c15 c16 c17 c18;
lb=[0 0 0 ];
ub=[10 10 20];
x0=[0 0 0];
H=[2 1 0
   -2 -1 2
   1 1 1
   -1 -1 1];
I=[(y+2*z+k)
   (0.5*x-y-2*z)
   m
   (0.5*x+3.76*z+k)];
options= optimset('MaxIter',5000,'MaxFunEvals',40000);
[XY,fval,exitflag,output] = fmincon('Objfcn',x0,H,I,[],[],lb,ub,'consfcn',options)


where 'Objfcn' is,

function [f]=Objfcn(x)
global c11 c12 c13 c14 c15 c16 c17 c18;
a=x(1);
b=x(2);
d=x(3);
f=sum([c11*a*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b^2 , c12*(c16 + 2*c17 + c18 - 2*a - b)*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a+b),c13*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b)^2 - d*(a + b + 1/2*c15 - d + 94/25*c17 + c18)].^2);


and 'confscn' is,

function [c,ceq] = consfcn(x)
global c11 c12 c13 c14 c15 c16 c17 c18;
a=x(1);
b=x(2);
d=x(3);
c=[];
ceq=[c11*a*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b^2 ; c12*(c16 + 2*c17 + c18 - 2*a - b)*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a+b);c13*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b)^2 - d*(a + b + 1/2*c15 - d + 94/25*c17 + c18)];


kindly do some comments. because this code gives good results only with certain values of constant. most of the time 'exit flag' gives '-2', and u know what does it mean :-)

takecare

"rabbahs " <shabbarraza2000@yahoo.com> wrote in message <hoik3g$qi0$1@fred.mathworks.com>...
> Walter Roberson <roberson@hushmail.com> wrote in message <ho7ejp$hro$1@canopus.cc.umanitoba.ca>...
> > rabbahs wrote:
> >
> > > can u give me the expression to calculate the positive value of a,b and d ?
> >
> > Well, as I wrote before:
> >
> > ====
> > The solution process involved reducing it down to a problem in 3
> > expressions, with the other 5 reducing to 0 by linear combinations of
> > the 3 remaining variables. Those expressions (which have to be solved to
> > equal 0) are:
> >
> > c11*a*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b^2,
> >
> > c12*(c16 + 2*c17 + c18 - 2*a - b)*(a + b + 1/2*c15 - d + 94/25*c17 +
> > c18) - b*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a+b),
> >
> > c13*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b)^2 - d*(a + b + 1/2*c15 - d +
> > 94/25*c17 + c18)
> >
> > As before, c18 is what you wrote as the constant c1(8) and so on.
> >
> > The solution to this trio of equations involves finding the roots of a
> > quartic.
> > =====
> >
> > So substitute your known constants into the three equations above, and
> > then proceed with traditional elimination. For example, d appears in the
> > first one multiplied only by a and a constant, so move the -b^2 to the
> > right hand side, divide by c11*a, and you then have a linear expression
> > on the left hand side that can be used to isolate d. That d value can
> > then be substituted into the two remaining equations. The second
> > expression will then become cubic in b but the third expression will
> > become quartic in b, so it would probably be easier to solve the second
> > expression first for a (make sure you consider all three roots), and
> > then use that to solve the third expression.
> >
> >
> > After substitution of the constants, the three expressions will be of
> > the general form
> >
> > E1*a*(a+b+E2-d)/E3-b^2,
> >
> > E4*(E5-2*a-b)*(a+b+E2-d)/E6-b*(E7-2*d+2*a+b),
> >
> > E8*(E9-2*d+2*a+b)/E10-d*(a+b+E2-d)
> >
> > where the E's are a combination of constants determined by the c1(:)
> > values. There is a general symbolic solution to this, but the general
> > solution has a repeated subexpression a quartic polynomial. In Matlab
> > notation, this quartic can be solved via
> >
> > P = roots([E10*(E1+4*E3) * (E1*E6*E4+E1*E6^2-E3*E4^2),
> >
> > E1*(8*E8*E4*E3*E6 - 2*E10*E2*E4^2*E3 - 4*E10*E3*E6*E4*E5 -
> > E1*E10*E6*E7*E4 - 4*E10*E6*E7*E4*E3 + E1*E8*E4^2 + 3*E1*E10*E2*E6*E4
> > + 4*E10*E2*E4*E3*E6 + 2*E1*E8*E4*E6 - E1*E10*E6*E4*E5 +
> > 2*E1*E10*E2*E6^2 + 4*E8*E4^2*E3),
> >
> > E1*(E1*E10*E6*E7*E4*E5 + 4*E1*E8*E2*E4*E6 - 2*E1*E8*E4*E5*E6 +
> > E10*E3*E4^2*E5^2 + 2*E1*E10*E2*E6^2*E7 + 4*E1*E8*E6^2*E7 -
> > E1*E10*E6^2*E7^2 + 2*E1*E8*E6*E7*E4 - 4*E8*E9*E1*E6^2 -
> > 4*E8*E4^2*E3*E5 + 2*E1*E10*E2^2*E6*E4 - E8*E9*E1*E4^2 +
> > 2*E1*E8*E2*E4^2 + 2*E10*E2*E4^2*E3*E5 - 3*E1*E10*E2*E6*E4*E5 -
> > 2*E1*E8*E4^2*E5 - E1*E10*E2*E6*E7*E4 - 4*E8*E9*E1*E4*E6),
> >
> > E1^2*E4*E5*(-E8*E4*E5 - 2*E4*E8*E9 + 4*E8*E2*E4 + 2*E8*E6*E7 -
> > 4*E8*E9*E6 + 2*E10*E2^2*E6 - E10*E2*E6*E7 + 4*E8*E2*E6),
> >
> > -E8*E9*E1^2*E4^2*E5^2 + 2*E8*E1^2*E2*E4^2*E5^2 )
> >
> > where the various E's are copied by inspection of the coefficients of
> > the three equations above after substitution of the constants.
> >
> > Once you have the roots of this P value, substitute them each in turn into:
> >
> > a = -P*E3*((-E4 - 2*E6)*P + E4*E5) /
> > ((E1*E6 - 2*E3*E4)*P - E6*E7*E1 + 2*E6*E1*E2)
> >
> > b = P
> >
> > d = (-(E1 + 4*E3)*(E1*E6*E4 + E1*E6^2 - E3*E4^2)*P^3 - E1*(3*E1*E4*E6*E2
> > - E1*E4*E5*E6 - 2*E2*E4^2*E3 - E1*E4*E6*E7 - 4*E3*E6*E4*E5 +
> > 4*E2*E6*E4*E3 + 2*E1*E6^2*E2 - 4*E6*E7*E4*E3)*P^2 - E1*(E3*E4^2*E5^2 -
> > E1*E6^2*E7^2 - 3*E1*E4*E5*E6*E2 + 2*E1*E2*E6^2*E7 + 2*E1*E2^2*E4*E6 -
> > E1*E2*E4*E6*E7 + E1*E4*E5*E6*E7 + 2*E2*E4^2*E5*E3)*P -
> > E1^2*E2*E4*E5*E6*E7 + 2*E1^2*E2^2*E4*E5*E6) /
> > (((-E4 - 2*E6)*P + E4*E5)*E1*((E1*E6 - 2*E3*E4)*P - E6*E7*E1 +
> > 2*E6*E1*E2)
> >
> > and see if you can find a root set of P that makes a, b, and d all
> > non-negative. If none of the root-sets work, then there is no hope with
> > that combination of constants. (This set of expressions was found via
> > the Maple symbolic combination engine -- I don't try to do these things
> > by hand!)
> >
> > Once you have an all-positive (a, b, d) tuple, substitute the respective
> > constants and a, b, d values into
> >
> > c = 1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b
> > e = 94/25*c17
> > f = c16 + 2*c17 + c18 - 2*a - b
> > g = c14 - a - b - d
> > p = a + b - d + c15/2 + 94/25*c17 + c18
> >
> > and see if you get all-positive results. If so, you are done. If not,
> > then go back and try the next root-set of the output of the roots()
> > call. If none of the four work, then you cannot solve the problem with
> > that combination of constants.
> >
> > I do not recommend trying to solve this system via a minimizer: the
> > constants either have the right relationship for the exact solution, or
> > they do not have the right relationship and no amount of minimization
> > would fix that. The only time I might suggest a minimizer on this is if
> > the "constants" are actually intervals and you are hoping that somewhere
> > in the intervals is a solution such that all the variables are positive.
> >
> >
> > If you explore this system of equations further, you will find that in
> > general it has singularities. The singularities _tend_ to be in the
> > negative range, but with some combinations of values, the singularities
> > can occur in the positive ranges. For example, with the set of random
> > c1(:) constants I indicated before that I posted with, except with c18
> > left undetermined, I found that the graph was quite wobbly between 8.1
> > and 8.4 with multiple roots _and_ a singularity within a range of about
> > 0.102. The singularity itself, Matlab cannot isolate, because a change
> > in the last bit causes a jump from about -21 to about +7: there might
> > actually be a zero in there, if you used used much higher precision...
> > but be off by .1 and the values jump between +/- 10^64 !
>
> Thanks Mr. Walter,
>
> currently i m using fmincon command to solve this problem. I reduced the equations in to 3 set of equations (with variable a,b and 9, as u advise earlier) and that try to bound my solution with the constrains available with this function.
>
> here are some details,
>
> global c11 c12 c13 c14 c15 c16 c17 c18;
> lb=[0 0 0 ];
> ub=[10 10 20];
> x0=[0 0 0];
> H=[2 1 0
> -2 -1 2
> 1 1 1
> -1 -1 1];
> I=[(y+2*z+k)
> (0.5*x-y-2*z)
> m
> (0.5*x+3.76*z+k)];
> options= optimset('MaxIter',5000,'MaxFunEvals',40000);
> [XY,fval,exitflag,output] = fmincon('Objfcn',x0,H,I,[],[],lb,ub,'consfcn',options)
>
>
> where 'Objfcn' is,
>
> function [f]=Objfcn(x)
> global c11 c12 c13 c14 c15 c16 c17 c18;
> a=x(1);
> b=x(2);
> d=x(3);
> f=sum([c11*a*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b^2 , c12*(c16 + 2*c17 + c18 - 2*a - b)*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a+b),c13*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b)^2 - d*(a + b + 1/2*c15 - d + 94/25*c17 + c18)].^2);
>
>
> and 'confscn' is,
>
> function [c,ceq] = consfcn(x)
> global c11 c12 c13 c14 c15 c16 c17 c18;
> a=x(1);
> b=x(2);
> d=x(3);
> c=[];
> ceq=[c11*a*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b^2 ; c12*(c16 + 2*c17 + c18 - 2*a - b)*(a + b + 1/2*c15 - d + 94/25*c17 + c18) - b*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a+b);c13*(1/2*c15 - 2*d - c16 - 2*c17 + 2*a + b)^2 - d*(a + b + 1/2*c15 - d + 94/25*c17 + c18)];
>
>
> kindly do some comments. because this code gives good results only with certain values of constant. most of the time 'exit flag' gives '-2', and u know what does it mean :-)
>
> takecare



********************************************************************

where,

c14=m;c15=x;
c16=y;c17=z;c18=k

hope it makes sense !

rabbahs wrote:

> Thanks Mr. Walter,
>
> currently i m using fmincon command to solve this problem. I reduced the
> equations in to 3 set of equations (with variable a,b and 9, as u advise
> earlier) and that try to bound my solution with the constrains available
> with this function.

> kindly do some comments. because this code gives good results only with
> certain values of constant. most of the time 'exit flag' gives '-2', and
> u know what does it mean :-)

It means what I said before: that I would not recommend trying to solve
the problem with a minimizer. I showed the form of the exact solution of
the problem, and if you have an exact solution, you do not need a
minimizer: you just plug in the numbers, and what comes out either meets
the constraints or it doesn't (implying that searching for a solution
would be useless.)

I could recompute the exact solution (sorry, I deleted it yesterday...
on the other hand I might well be able to recover it from my Time
Machine backups) and email it to you in Maple format, but it was about
12 megabytes. Or I could convert it to an optimized Matlab procedure (or
at least could do so one variable at a time), but that form is even
longer and is totally incomprehensible.

The only point in using a minimizer would be if some of your constants
are not firm values and you are looking for the least change in your
"constants" that would allow the problem to be solved. To take this
approach, you would have to define a penalty function that trades off
the adherence to the original value of the constants against the work
required to find a solution, and which might perhaps weight the
different "constants" differently... different "stiffnesses" for each.

Walter Roberson <roberson@hushmail.com> wrote in message <hoinjr$r39$1@canopus.cc.umanitoba.ca>...
> rabbahs wrote:
>
> > Thanks Mr. Walter,
> >
> > currently i m using fmincon command to solve this problem. I reduced the
> > equations in to 3 set of equations (with variable a,b and 9, as u advise
> > earlier) and that try to bound my solution with the constrains available
> > with this function.
>
> > kindly do some comments. because this code gives good results only with
> > certain values of constant. most of the time 'exit flag' gives '-2', and
> > u know what does it mean :-)
>
> It means what I said before: that I would not recommend trying to solve
> the problem with a minimizer. I showed the form of the exact solution of
> the problem, and if you have an exact solution, you do not need a
> minimizer: you just plug in the numbers, and what comes out either meets
> the constraints or it doesn't (implying that searching for a solution
> would be useless.)
>
> I could recompute the exact solution (sorry, I deleted it yesterday...
> on the other hand I might well be able to recover it from my Time
> Machine backups) and email it to you in Maple format, but it was about
> 12 megabytes. Or I could convert it to an optimized Matlab procedure (or
> at least could do so one variable at a time), but that form is even
> longer and is totally incomprehensible.
>
> The only point in using a minimizer would be if some of your constants
> are not firm values and you are looking for the least change in your
> "constants" that would allow the problem to be solved. To take this
> approach, you would have to define a penalty function that trades off
> the adherence to the original value of the constants against the work
> required to find a solution, and which might perhaps weight the
> different "constants" differently... different "stiffnesses" for each.

********************************************************************

i am agree with your point that if we have exact solution to the problem then v dont have to use the minimization function.
If u recover the solution then kindly let me know

tc

rabbahs wrote:
> Walter Roberson <roberson@hushmail.com> wrote in message
> <hoinjr$r39$1@canopus.cc.umanitoba.ca>...
>> rabbahs wrote:

>> It means what I said before: that I would not recommend trying to
>> solve the problem with a minimizer. I showed the form of the exact
>> solution of the problem, and if you have an exact solution, you do not
>> need a minimizer: you just plug in the numbers,

>> I could recompute the exact solution (sorry, I deleted it yesterday...
>> on the other hand I might well be able to recover it from my Time
>> Machine backups) and email it to you in Maple format, but it was about
>> 12 megabytes.

Actually it was 19 megabytes of incomprehensible masses of symbols.

> If u recover the
> solution then kindly let me know

I will not send the final solution because it is just too big. Instead, I will
  show how to get the solution in Maple. I highly recommend using the
command-line version of Maple for this purpose.


First, reduce down to the 3 equations in a, b, d, like we discussed before. I
will refer to this as eqmat3.

eqmat3sol := solve(eqmat3,[a,b,d]) assuming real;

OR//

eqmat_sol := RealDomain:-solve(eqmat3,[a,b,d]);

If you are using the command line version of Maple, it will print out a single
solution of the form

[[a = mess#1 involving %1, b = %1, d = mess#2 involving %1]]

and then it will print out

%1 := mess#3

This means that Maple has identified a significant common sub-expression to
several terms, and has broken it out so as to make the output easier to read.
If you were using one of the GUI versions of Maple, you would have to go
through a bit of trouble to get Maple to print things out with this common
sub-expression made clear. The common sub-expression will be a RootOf() call
close to 60 lines long that sets out a quartic equation. There will be four
roots to the quartic, and they will all be terribly long and messy; you will
have a fair bit of trouble understanding the solution if you proceed to ask
Maple to find all four roots of the quartic. And it will take a long time,
too. So we will step around the issue a bit to make things faster.

Now, command:

P1val := %1:

This isolates the RootOf() into a variable. Now construct a variant of the
solution matrix rewritten with a symbolic variable for this RootOf:

eqmat3_solP1 := subs(P1val=P1,eqmat3_sol):

Once we have found candidate solutions for the RootOf(), we will go back and
substitute them in individually for P1, and ensure that a and d come out
non-negative (b is the same as the root, and since b must non-negative as
well, we will only select non-negative roots.) If a and d pass, then we go
back and substitute the a, b, and d into the other 5 variables that have been
rewritten in terms of a, b, and d, as shown earlier in the discussion.

The next step towards finding the solution is formally necessary according to
the Maple documentation, but I've never seen a problem if it isn't used. It
certainly doesn't hurt and probably makes things faster, so we might as well
do it anyhow:

Roc := collect(op(P1val),_Z):

C4v := factor(coeff(Roc,_Z,4)):
C3v := factor(coeff(Roc,_Z,3)):
C2v := factor(coeff(Roc,_Z,2)):
C1v := factor(coeff(Roc,_Z,1)):
C0v := factor(coeff(Roc,_Z,0)):

This isolates the coefficients of the _Z terms of the RootOf.

Now the time saving step: instead of finding the RootOf the big messy
equation, we will find the RootOf a generic quartic, and substitute in the
above coefficients afterwards.

G4 := RootOf(C4*_Z^4+C3*_Z^3+C2*_Z^2+C1*_Z+C0):
G4S := [allvalues(G4)]:

I pre-computed G4S and saved it away so that I don't have find the solutions
over and over.

tmpproc := unapply(G4S[1],[C0,C1,C2,C3,C4]):
G4S1proc := codegen[optimize](tmpproc):
tmpproc := unapply(G4S[2],[C0,C1,C2,C3,C4]):
G4S2proc := codegen[optimize](tmpproc):
tmpproc := unapply(G4S[3],[C0,C1,C2,C3,C4]):
G4S3proc := codegen[optimize](tmpproc):
tmpproc := unapply(G4S[4],[C0,C1,C2,C3,C4]):
G4S4proc := codegen[optimize](tmpproc):

After this set of Maple commands, G4S*proc are procedures for efficiently
computing the individual quartic roots given the coefficients. Again, you may
wish to save these procedures away for future use.

Now:

P1sol1 := G4S1proc(C0,C1v,C2v,C3v,C4v):
P1sol2 := G4S2proc(C0,C1v,C2v,C3v,C4v):
P1sol3 := G4S3proc(C0,C1v,C2v,C3v,C4v):
P1sol4 := G4S4proc(C0,C1v,C2v,C3v,C4v):

These are now the individual roots for P1, the RootOf(), which is also the
value for b.

If you want to proceed symbolically from here, you could construct entire
solution sets using (e.g.)

eval(eqmat3_solP1, P=P1sol1);

but you are going to get messy.

I would recommend this as being about the point that you switch from symbolic
to numeric. Code as far as the P1sol* in terms of the constants c1(1), c1(2)
and so on (which I refer to as c11, c12 and so on here), and then code the
eqmat3_solP1 in terms of the constants and P1, and then code the rest of the
variables in terms of the a, b, and d that you get out. Only bother to plop in
the values of P1 that came out nonnegative and real; only bother to continue
on to the extended equations if a and d also come out nonnegative and real;
and, of course, afterwards sanity-check everything to be sure that all the
variables came out non-negative and real.

With this construction, if any of the variables come out negative or complex,
then there is no solution for that set of c1 constants.

If you are going to be using a number of different c1 constant sets, then you
might want to take the extra step of:

tmpproc := unapply(P1sol1,[c11, c12, c13, c14, c15, c16, c17, c18]):
P1sol1proc := codegen[optimize](tmpproc):
tmpproc := unapply(P1sol2,[c11, c12, c13, c14, c15, c16, c17, c18]):
P1sol2proc := codegen[optimize](tmpproc):
tmpproc := unapply(P1sol3,[c11, c12, c13, c14, c15, c16, c17, c18]):
P1sol3proc := codegen[optimize](tmpproc):
tmpproc := unapply(P1sol4,[c11, c12, c13, c14, c15, c16, c17, c18]):
P1sol4proc := codegen[optimize](tmpproc):

(The roots should come out pretty similar, actually -- just like a quadratic
has two solutions because it has a basic formula with one place to choose +/-,
a quartic has four solutions because it has a basic formula with two places to
choose +/- .)