Solving 4 non-linear equations for 4 unknowns using fsolve

Hello,

I am trying to fit 2 non-linear equations (dashed, blue lines) to 2 sets of experimental data (solid, red or black lines). Within these 2 equations are 4 unknowns (u1, n1, u2, and n2). I have therefore developed 4 equations such that the endpoints of my new equations should be equal those of the experimental data. Now I have 4 equations and 4 unknowns. Here are the issues that I am experiencing using fsolve:

1. When plugging in the 4 solutions back into my 2 non-linear equations, the endpoints do not match the experimental data. I expect that the endpoint values should be the same because I set these in fsolve_function.
2. When making the initial guesses for the solutions, 2 of the solutions (e.g. n1 and n2) are always the same as my initial guess, and the other 2 solutions (e.g. u1 and u2) change. I expect all 4 values to change to values that are different from my initial guesses.

What could be the causes?

Here is my function:

function F = fsolve_function (Y)

%Unpack variables
u1 = Y(1);
n1 = Y(2);
u2 = Y(3);
n2 = Y(4);

%Constants
e = 1.602 *10^-19;  %Coulomb, electron charge

%Constraints from experimental data
A = 0.000438;       %value at B=0 for Gxx
B = 0.0003;         %value at B=14 for Gxx
C = 0.0000061841;   %slope at B=0 for Gxy
D = 0.0000381;      %value at B=14 for Gxy


%Input equations

%Gxx, inserting value at B=0T
equation1 = n1*e*u1 + n2*e*u2 - A; 

%Gxx, inserting value at B=14T
equation2 = ( (n1*e*u1) / (1+u1^2*(14^2)) ) + ( (n2*e*u2) / (1+u2^2*(14^2)) ) - B; 

%Gxy, inserting slope at B=0T
equation3 = n1*e*u1^2 + n2*e*u2^2 - C; 

%Gxy, inserting value at B=14T
equation4 = ( (n1*e*u1^2*14) / (1+u1^2*14) ) + ( (n2*e*u2^2*14) / (1+u2^2*14) ) - D; 

%Pack up our function
F = [equation1; equation2; equation3; equation4; ];

Here is my script:

%Initial guesses
u1o = 0.009; 
n1o = 5 * 10^17; 
u2o = 0.5; 
n2o = 2 * 10^19; 
Y0 = [u1o; n1o; u2o; n2o];

%Call fsolve function
fhandle = @fsolve_function;
Y = fsolve (fhandle, Y0);

%Display solutions for our 4 variables
disp ('u1, n1, u2, n2');
disp (Y(1:4));

Here are my results: As one can see. The endpoints of the simulated data (dashed, blue lines) do not match those of the experimental data (solid red or black lines) even though I set these values in my fsolve_function. Secondly the values for n1 and n2 are the same as my initial gueses while the others (u1 and u2) changes. If someone could help me understand perhaps what the issue may be that would be very helpful!

Thanks,

Chris

[Walter Roberson]

There are two real-valued solutions to the equations:

Y1 = 0.1437438700474650, Y2 = 4.82587963631918*10^15, Y3 = -0.0299508173739704, Y4 = -6.812474451761340*10^16

Y1 = -0.02995081737398195, Y2 = -6.812474451761304*10^16, Y3 = 0.1437438700474649, Y4 = 4.825879636317568*10^15

The other 6 solutions all involve complex values.

Answer 1

[Matt J]

The modified code below yields a numerical solution decently close to the second analytical solution given by Walter. The basic problem was that you measured n1 and n2 in terms of numbers of electrons. This is not a good choice of units for the problem, because the equation function is very insensitive to increases/decreases by 1 electron as compared to unit changes in the other parameters, u1 and u2. By instead expressing the equations in terms of total charge n1e and n2e in Coulombs, the sensitivity to the different parameters balances out much better.

%Initial guesses
u1o = 0.009; 
n1o = 5e-2; 
u2o = 0.5; 
n2o = 2; 
Y0 = [u1o; n1o; u2o; n2o];

%Call fsolve function
fhandle = @fsolve_function;
opts=optimoptions('fsolve','Display','iter','TolFun',1e-12);
[Y,fval,exitflag] = fsolve (fhandle, Y0,opts);

%Constants
e = 1.602 *10^-19;  %Coulomb, electron charge

Y([2,4])=Y([2,4])/e;

fval,exitflag

%Display solutions for our 4 variables
disp ('u1, n1, u2, n2');
C=(num2cell(Y));
C{:}


function F = fsolve_function (Y)

    %Unpack variables
    u1 = Y(1);
    n1e = Y(2);
    u2 = Y(3);
    n2e = Y(4);


    %Constraints from experimental data
    A = 0.000438;       %value at B=0 for Gxx
    B = 0.0003;         %value at B=14 for Gxx
    C = 0.0000061841;   %slope at B=0 for Gxy
    D = 0.0000381;      %value at B=14 for Gxy


    %Input equations

    %Gxx, inserting value at B=0T
    equation1 = n1e*u1 + n2e*u2 - A; 

    %Gxx, inserting value at B=14T
    equation2 = ( (n1e*u1) / (1+u1^2*(14^2)) ) + ( (n2e*u2) / (1+u2^2*(14^2)) ) - B; 

    %Gxy, inserting slope at B=0T
    equation3 = n1e*u1^2 + n2e*u2^2 - C; 

    %Gxy, inserting value at B=14T
    equation4 = ( (n1e*u1^2*14) / (1+u1^2*14) ) + ( (n2e*u2^2*14) / (1+u2^2*14) ) - D; 

    %Pack up our function
    F = [equation1; equation2; equation3; equation4; ];

end

This gives me:

u1, n1, u2, n2

ans =

   -0.0300


ans =

  -6.8125e+16


ans =

    0.1437


ans =

   4.8259e+15

[Walter Roberson]

Analytically, Y1 is the set of roots of a degree 8 polynomial, and the other variables are degree 7 polynomials in Y1.

... not nice polynomials, but polynomials. For example.

R = Y1;
Y2 = 8755012204935750201112908752536978151288890177813486758257767702284373502848255044532791108125066/37094697431329030102201238269213882491399313262077363526626133919642866755656767620183169693 * R^7 - ...
  378856539749416824412013428047210835599190735215388865104650782583751582785768716444349145953411/37094697431329030102201238269213882491399313262077363526626133919642866755656767620183169693     * R^6 + ...
  98432853485230761182005542408267981001337433399456227311391463628457116485894713846155150720960/37094697431329030102201238269213882491399313262077363526626133919642866755656767620183169693      * R^5 - ...
  27962962520599252583762104859734641779889928285945889630711273923541165197225318905221906320880/37094697431329030102201238269213882491399313262077363526626133919642866755656767620183169693      * R^4 - ...
  308216035260224196829655287907072958412694771400759109927852662446634904950119443094400644046/37094697431329030102201238269213882491399313262077363526626133919642866755656767620183169693        * R^3 - ...
  108441310057192148021222082274451632121192927877953892372277908362726479517268511272042786860/37094697431329030102201238269213882491399313262077363526626133919642866755656767620183169693        * R^2 + ...
  8350417447337202706363180983130138572802372532652765840049571119574384765702151026420948715/37094697431329030102201238269213882491399313262077363526626133919642866755656767620183169693          * R   + ...
  111728589093862187739186415682620972846352251032216174804370977345691174076457888917331885/37094697431329030102201238269213882491399313262077363526626133919642866755656767620183169693

[Walter Roberson]

Oh yes, and Y1 is the set of roots of

68270502110574346712672641163215073783648026173556493777530875649523270689660862464 * Z^8 - ...
2855529907609170739195229222548832307587949808189604607194299298563103294713430016  * Z^7 + ...
784440707766226545188364555997666186250357690273799349031332505557438893254091232   * Z^6 - ...
216466632607974965268412283559078238219659271767095852323723699764179685008756736   * Z^5 - ...
2540241169862984402955650135774412878075901682320618381658630012661512217359116     * Z^4 - ...
882316986721974487685827460065764355079676423033507749733197764209534536162304      * Z^3 + ...
57843594902826116320965552582878306746006826402495469547105052988089617647748       * Z^2 + ...
753942968284174553979797778926098357020195931184331774827163392982334665216         * Z   - ...
56858285817782871654260116722799255821670781443196834071920926362838424443

"Obviously".

[Christian Greenhill]

Thank you Matt and Walter. Your modified script works!