Fsolve-positive answers only?

Hi, I am a new Matlab user so I am very inexperienced with the program...yet. Hope I will get a hang of it soon.

I have non-linear multivariable equations:

function F =root8d(x)
F(1)=0.02-x(1)-x(4)-x(5)-x(6)-x(7)-2*x(8);
F(2)=10^-6-x(2)-x(4)-x(5)-x(6)-x(7)-2*x(8);
F(3)=0.01-x(3)-x(4)-2*x(5)-3*x(6)-4*x(7)-4*x(8);
F(4)=x(4)-10^8.35*x(1)*x(2)*x(3);
F(5)=x(5)-10^15*x(1)*x(2)*x(3)^2;
F(6)=x(6)-10^19.62*x(1)*x(2)*x(3)^3;
F(7)=x(7)-10^21.12*x(1)*x(2)*x(3)^4;
F(8)=x(8)-10^31.02*x(1)^2*x(2)^2*x(3)^4;

The problem is, when I use normal command

fun = @root8d;
x0 = [whatever number you put];
x = fsolve(fun,x0).

I DO get answers, but the problem is, these equations are from solution equilibrium so the answers has to be positive(or nonnegative at least). I know fsolve doesn't allow constraints.

Are there any ways to make the fsolve to yield positive answers only?

Thanks

Answer 1

[Walter Roberson]

fun2 = @(x) sum(root8d(x).^2);
A = [];
b = [];
Aeq = [];
beq = [];
lb = zeros(1,8);
ub = inf(1,8);
nonlcon = [];
options = ... you will probably need to increase iterations
one_solution = fmincon(fun2, x0, A, b, Aeq, beq, lb, ub, nonlcon, options);

[Walter Roberson]

solving for the first three variables analytically is easy. The fourth variable is a cubic, which gives you a long explicit expression. Everything goes down hill from there, with the fifth variable not practically tractable.

[Walter Roberson]

Ah, Maple says solving the first 5 simultaneously only gets you to a quintic. But beyond that, you have trouble because quintics have no analytic solution.

[Paul You]

Thank you Mr. Robertson. I will try these methods.

A graduate student working in same research group suggested me to use loop until I get positive numbers. What are the basic codes of setting up loops?

[Walter Roberson]

I would not recommend a loop for this purpose. fsolve() can end up examining the values in any order. The bounded fmincon() will only search the positive area.

[Walter Roberson]

My tests indicate that there is a solution very close to [0.019998999999661933, 9.10591091655229996e-19, 0.00999676207120531179, 5.92705032715569534e-13, 1.81987230951624782e-09, 7.58419642443709219e-07, 2.39755731650481807e-07, 1.98656134362228394e-12]

The residue at that location is on the order of 6E-24 which is effectively down in the numeric noise range (that is, could vary notably if the operations were done in a different order.)

Knowing a good starting point is important to finding a solution. The starting point [0.02, 0, 0.01, 0, 0, 0, 0, 0] is pretty good.

[Walter Roberson]

You will find that one of the important steps in efficient searching is to vectorize your code. That can be done by using

function F =root8d(x1, x2, x3, x4, x5, x6, x7, x8)
F1 = 0.02 - x1 - x4 - x5 - x6 - x7 - 2 .* x8;
F2 = 10^(-6) - x2 - x4 - x5 - x6 - x7 - 2 .* x8;
F3 = 0.01 - x3 - x4 - 2*x5-3 .* x6 - 4 .* x7 - 4 .* x8;
F4 = x4 - 10^8.35 .* x1 .* x2 .* x3;
F5 = x5 - 10^15 .* x1 .* x2 .* x3.^2;
F6 = x6 - 10^19.62 .* x1 .* x2 .* x3.^3;
F7 = x7 - 10^21.12 .* x1 .* x2 .* x3.^4;
F8 = x8 - 10^31.02 .* x1.^2 .* x2.^2 .* x3.^4;
sh = [1, size(x1)];
F = [reshape(F1,sh); reshape(F2,sh); reshape(F3,sh); reshape(F4,sh); reshape(F5,sh); reshape(F6,sh); reshape(F7,sh); reshape(F8,sh)];

together with

residue_s = @(x1,x2,x3,x4,x5,x6,x7,x8) reshape(sum(root8d(x1,x2,x3,x4,x5,x6,x7,x8).^2, 1), size(x1));
residue_v = @(x) residue_s(x(1), x(2), x(3), x(4), x(5), x(6), x(7), x(8));

The modified version of your code accepts inputs of arbitrary dimension and runs the calculations on all of the items simultaneously, and splices the 8 pieces together along the first dimension (so if you passed in 93 x 17 matrices, you would get out 8 x 93 x 17). residue_s runs the function and converts those values into sum-of-squares residues the same shape as the original value. Together, these allow you to take approaches such as using ndgrid() to construct your input matrices and to get out one residue result for each location. This allows you to sweep over a parameter range looking for good minima for further exploration.

The residue_v function accepts a vector of 8 values and runs the residue calculation by passing them in to the appropriate slots in the vectorized calculation. This residue_v version can be used for situations such as fmincon or fminsearch that need to pass in vectors of parameters.

I find that fmincon with default parameters does not necessarily do well with this situation -- it can end up giving outputs higher than its initial point :(

[Paul You]

Thanks so much Mr Roberson..what are: A = []; b = []; Aeq = []; beq = []; lb = zeros(1,8); ub = inf(1,8); nonlcon = [];

[Walter Roberson]

fmincon() takes a series of positional parameters representing constraints. The plain upper and lower bound constraints (e.g., to say that a value needs to be non-negative, which is to say must be in the range 0 to +inf) happen to be defined to occur later on in the sequence.

The first pair of constraints that fmincon allows are the linear inequalities, of the form A*x <= b where A is a constant matrix and b is a constant vector. If you do not have linear inequalities, pass empty matrices instead.

The second pair of constraints that fmincon allows are the linear equalities, of the form Aeq*x == beq where A is a constant matrix and b is a constant vector. If you do not have linear equalities then pass empty matrices instead.

The third pair of constraints that fmincon allows are the lower bounds (lb) and upper bounds (ub), each constant vectors. If you do not have lower bounds and upper bounds, then pass empty matrices instead.

After that, fmincon allows a function handle that defines the nonlinear constraints. If you do not have nonlinear constraints then pass an empty matrix instead.

The order is A, b, Aeq, beq, lb, ub, nonlcon .

Instead of just passing in the right number of [] such as

fmincon(fun2, x0, [], [], [], lb, ub, [])

I prefer to assign those [] to variables using the usual names, so that anyone else reading the code can immediately see what I am passing in and does not need to count the [] to figure out what the parameters are intended to mean.

Answer 2

[Matt J]

FSOLVE doesn't allow bounds, but LSQNONLIN does

x = lsqnonlin(fun,x0,zeros(size(x0)));