MATLAB: The result of fsolve is wrong

fsolvesolve

I have system of linear equations in three unknowns.
  1. x(1)-1+(1-x(3))^(M1-1)=0
  2. x(2)-1+(1-x(3))^(M2-1)=0
  3. x(3)-(1+x(1)^(R1+1)+I+J)*e*Pi=0
However, the results in my Matlab is wrong.
How can I do to optimize this code? Or what's wrong in my codes?
Following is my Matlab code:
function T=Analysis0828(R1,R2,W1,W2,M,a1,a2,e)
x = sym('x',[1 2 3]);
%Equations  
M1=(M/(a1/a1+a2))
M2=(M/(a2/a1+a2))
G=((W1+1)/2)*(x(1)*(1-x(1)^(R1)))/(1-x(1));
H=(x(1)^(R1+1))*((W2+1)/2)*x(2)*(1-x(2)^(R2))/(1-x(2));
Pi=(1+e+(x(1)^(R1+1))*e+(G+H)*e)^(-1);
I=x(1)*(1-x(1)^(R1))/(1-x(1));
J=x(1)*(1-x(2)^(R2))/(1-x(2));
my_eqns = [x(1)-1+(1-x(3))^(M1-1); 
           x(2)-1+(1-x(3))^(M2-1);  
           x(3)-(1+x(1)^(R1+1)+I+J)*e*Pi];  
  g = matlabFunction(my_eqns);
    F = my_vectorize(g);
x0 = [0.4,0.5,0.01];
s = fsolve(F,x0);
T=(M1*M1*(1-s(1))*s(3)+M2*M2*(1-s(2))*s(3))/M
my_vectorize funtion code is as follows
 function f_vec = my_vectorize(fun)
    f_vec = @fun_vectorize;
    function out = fun_vectorize(x)
        x   = num2cell(x);
        out = fun(x{:});
    end
 end
When inputs is (3, 5, 5, 2, 10, 1, 1, 1)
And my results are
s =
    0.9557    0.5413    0.9557

Best Answer

We do not know what your my_vectorize does; in particular we do not know how it turns the function of three inputs into a function of one input. The order you arrange the inputs is important, in order for us be able to calculate which of the 37 complex-valued solutions is "closest" to [0.4,0.5,0.01]
Have you considered using the 'vars' option to matlabFunction instead of using your own custom vectorization function?
If we assume that the order for [0.4,0.5,0.01] is x1_1, x2_1, x1_2 (the first three entries of the 1 x 2 x 3 matrix of x symbols) then the one of the 37 solutions that has the smallest norm to [0.4, 0.5, 0.01] is
However, your use of the initial conditions [0.4,0.5,0.01] would imply that you want the solution that is "closest" to [0.4,0.5,0.01], and the closest is
x1_1 = 0.95573206349551101631165448445161
x2_1 = 0.95573206349551101631165448445161
x1_2 = 0.54130678130277676362675222900352
which turns out to be a re-arrangement of what you said you got as a solution. Perhaps your vectorization uses a different order.
When given the wrong order of initial conditions, fsolve() will wander off into complex values; when given the correct order of initial conditions that matches the variable in the x1_2 position being 0.01, then fsolve() will produce the numeric values 0.955732064923563 0.955732064923563 0.541306780698535
So I think that fsolve() is working -- but you have to be more careful about how you vectorize and about the order of initial conditions.
Related Question