How can I solve this system? fsolve solver fails

Hi, I have a system of equations that I need to solve, I've been using the fsolve function with a variety of initial guesses and increased the MaxFunEvals and MaxIter to no avail. Is there anything else I can try?

My system is given by:

function eq3 = equilibrium3(x,S,P,A,B)
a = x(1); b = x(2); c = x(3); d = x(4); e = x(5); f = x(6);
Sx=S(1);Sy=S(2);Sz=S(3);
i1=ICM([Sx+min(Sx,Sy)+1,max(0,Sy-Sx),Sz-1],P); %scf x wins%
i2=ICM([max(0,Sx-Sy),Sy+min(Sx,Sy)+1,Sz-1],P); %scf y wins%
i3=ICM([Sx+1.5,Sy-0.5,Sz-1],P); %sff x wins%
if (Sy<=Sx)&&(Sz<=Sx)
    if Sy>Sz
        i4=P;
    else
        if Sz>Sy
            i4=[P(1) P(3) P(2)];
        else
            i4=[P(1) 0.5*(P(2)+P(3))*[1 1]];
        end
    end
else    
i4=ICM([Sx+min(Sx,Sy)+min(Sx,Sz),max(0,Sy-Sx),max(0,Sz-Sx)],P); %scoc x wins%
end
if (Sx<=Sy)&&(Sz<=Sy)
    if Sx>Sz
        i5=P;
    else
        if Sz>Sx
            i5=[P(3) P(1) P(2)];
        else
            i5=[0.5*(P(2)+P(3)) P(1) 0.5*(P(2)+P(3))];
        end
    end
else    
i5=ICM([max(0,Sx-Sy),Sy+min(Sx,Sy)+min(Sy,Sz),max(0, Sz-Sy)],P); %scoc y wins%
end
if (Sy<=Sz)&&(Sx<=Sz)
    if Sy>Sx
        i6=[P(3) P(2) P(1)];
    else
        if Sx>Sy
            i6=[P(2) P(3) P(1)];
        else
            i6=[0.5*(P(2)+P(3))*[1 1] P(1)];
        end
    end
else    
i6=ICM([max(0,Sx-Sz),max(0,Sy-Sz),Sz+min(Sx,Sz)+min(Sy,Sz)],P); %scoc z wins%
end
i7=ICM([Sx+min(Sx,Sz)+0.5,Sy-0.5,max(0,Sz-Sx)],P); %sfc x wins%
i8=ICM([max(0,Sx-Sz),Sy-0.5,Sz+min(Sx,Sz)+0.5],P); %sfc z wins%
i9=ICM([Sx,Sy+min(Sy,Sz),max(0,Sz-Sy)],P); %fsc y wins%
i10=ICM([Sx,max(0,Sy-Sz),Sz+min(Sy,Sz)],P); %fsc z wins%
i11=ICM([Sx,Sy+1,Sz-1],P); %fsf y wins%
i12=ICM([Sx,Sy-0.5,Sz+0.5],P);%ff z wins%
f1=@(hy) (i1(1).*H(a,A,B)+i2(1).*H(hy,A,B))./(H(a,A,B)+H(hy,A,B));
f2=@(hz) (i7(1).*H(a,A,B)+i8(1).*H(hz,A,B))./(H(a,A,B)+H(hz,A,B));
f3=@(hy,hz) (i4(1).*H(a,A,B)+i5(1).*H(hy,A,B)+i6(1).*H(hz,A,B))./(H(a,A,B)+H(hy,A,B)+H(hz,A,B));
f4=@(hx) (i2(2).*H(b,A,B)+i1(2).*H(hx,A,B))./(H(b,A,B)+H(hx,A,B));
f5=@(hx,hz) (i5(2).*H(b,A,B)+i4(2).*H(hx,A,B)+i6(2).*H(hz,A,B))./(H(b,A,B)+H(hx,A,B)+H(hz,A,B));
f6=@(hz) (i9(2).*H(c,A,B)+i10(2).*H(hz,A,B))./(H(c,A,B)+H(hz,A,B));
f7=@(hx,hy) (i6(3).*H(d,A,B)+i4(3).*H(hx,A,B)+i5(3).*H(hy,A,B))./(H(d,A,B)+H(hx,A,B)+H(hy,A,B));
f8=@(hx) (i8(3).*H(e,A,B)+i7(3).*H(hx,A,B))./(H(e,A,B)+H(hx,A,B));
f9=@(hy) (i10(3).*H(f,A,B)+i9(3).*H(hy,A,B))./(H(f,A,B)+H(hy,A,B));
eq3(1)=i3(1).*(1-b).*(1-e)+(1-d).*integral(f1,0,b)+(1-b).*integral(f2,0,e)+integral2(f3,0,b,0,d); %dE[X]/da%
eq3(2)=i7(2).*a+(1-d).*integral(f4,0,a)+integral2(f5,0,a,0,d); %dE[Y]/db%
eq3(3)=i12(2).*(1-a)+i11(2).*(1-a).*(1-f)+(1-a).*integral(f6,0,f); %dE[Y]/dc%
eq3(4)=a.*b+integral2(f7,0,a,0,b); %dE[Z]/dd%
eq3(5)=a.*(1-b)+(1-b).*integral(f8,0,a); %dE[Z]/de%
eq3(6)=c.*(1-a)+(1-a).*integral(f9,0,c); %dE[Z]/df%

[Walter Roberson]

Please give some sample S, P, A, B .

Your code appears to use a function H that you have not shown the code for.

Could you confirm that you are looking for a vector of 6 values, x, such that ideally eq3(1:6) become 0?

Are there any restrictions on the range of any of the inputs? For example can we assume they will all be re-valued and that you are expecting all of the intermediate calculations to be real valued?

Can you post the equations?

If you are trying to vary x with the other parameters given, then everything you have assigned from Sx up to i12 appears to be independent of your input values, so it would be a good idea from an efficiently point of view to pre-compute those and pass them in to the function.

Answer 1

[John D'Errico]

Lots of reasons why fsolve might fail.

1. Your functions appear to be rather severely non-differentiable. In fact, it may well be they have discontinuities. But fsolve requires differentiability, and certainly continuity.

2. Many systems of equations need not have an exact solution. Do you KNOW there is a solution? How do you know that as a fact? If not, then the lack of a solution by fsolve may be merely a reflection of the lack of a solution.

3. Cranking down on the tolerances, increasing maxfunevals, etc., these are silly things, but ones that people often try first. They are equivalent to stamping your feet, and demanding that fsolve try harder. Find a solution because I REALLY want it solved. All that happens is fsolve may take a little while longer to not find a solution. Instead, you would almost always be far better off trying different (hopefully better) sets of starting values.

4. Why do you think that fsolve has failed? Too often I've seen people worried when they see a message that the optimization has terminated. Of course, terminated is merely a synonym for done.

Answer 2

[halleyhit]

Solving f(x)=0 means the optimization of min(abs(f(x))). So, optimtool, which is a powerful toolbox in matlab may also help you.