How do I solve this system of non-linear equations using fsolve?

Good afternoon everyone,

I am facing an issue in my code where vpasolve takes too much time to solve a system of non-linear equations (2 equations and 2 unknowns). I want to use fsolve instead as it is much faster, however, everytime I change the initial guess in fsolve, the results change. This could be because the tolerances in fsolve are not tight enough, but I am unsure which tolerances / options I should be changing. The main part of the code is shown below (the constants are found in the attached file):

% System of Non-linear equations: (Trying to caluclate er1(1) and er1(2),
% the rest terms are constants)

f1 = @(er1) [u1^2 + w1^2 - ((er1(2)*c1_rb)^2 + (er1(1)*c1_jr)^2 - ...
        2*(er1(2)*c1_rb)*(er1(1)*c1_jr)*cos(pi + atan(pi*((1-...
        (er1(1)^2))^0.5)/(4*er1(1)) - atan(pi*((1-(er1(2)^2))^0.5)/...
        (4*er1(2))))));(((R1_ro/R1_j)*((L1_rb/L1_jr)^3)*((c1_jr/c1_rb)^2)*(((1 - er1(1)^2)/...
        (1 - er1(2)^2))^2)*(er1(2)/er1(1))*((sqrt((16 - pi^2)*...
        (er1(2)^2) + pi^2))/(sqrt((16 - pi^2)*(er1(1)^2) + pi^2))))-...
        (R1_ri/R1_j)) - ((((pi*(L1_rb^3)*R1_ro*(er1(2)^2))*(sqrt(((er1(2)^2)*...
        ((16/(pi^2))-1))+1))*(sin(atan(pi*((1-(er1(2)^2))^0.5)/...
        (4*er1(2)))))/(8*c1_rb*((1-(er1(2)^2))^2)))+(((R1_ro^3)*...
        L1_rb*2*pi)/(c1_rb*sqrt(1 - er1(2)^2)))-(pi*(R1_ri^2)*...
        (L1_jr^3)*(er1(1)^2)*(sqrt(((er1(1)^2)*((16/(pi^2))-1))+1))*...
        (sin(atan(pi*((1-(er1(1)^2))^0.5)/(4*er1(1)))))/(8*R1_j*c1_jr*...
        ((1-(er1(1)^2))^2)))+((R1_j*(R1_ri^2)*L1_jr*2*pi)/...
        (c1_jr*sqrt(1 - er1(1)^2))))/((pi*R1_ri*(L1_jr^3)*(er1(1)^2)*...
        (sqrt(((er1(1)^2)*((16/(pi^2))-1))+1))*(sin(atan(pi*((1-...
        (er1(1)^2))^0.5)/(4*er1(1)))))/(8*c1_jr*((1-(er1(1)^2))^2))) +...
        (R1_ri*(R1_j^2)*L1_jr*2*pi)/(c1_jr*sqrt(1 - er1(1)^2))))];

options = optimoptions('fsolve','Display','none');
    
z1 = fsolve(f1,[1e-8 1e-8],options);

er1_jr = z1(1);

er1_rb = z1(2);

I know roughly what the solution of the simultanous equations should be as I plotted (both equations for the range of er1(1) and er1(2) that I want, which is between 0 to 1) as shown in the image below:

The results I get do not match the values shown in the figure above, add to that the fact that the solution changes everytime I change the initial guess, as I mentioned above.

Let me know if you want me to clarify any part of my question.

Many thanks in advance!

[Zeyad]

I apologize for that, thanks for pointing it out. The file has been updated. R1_j = 5.5e-3

[Star Strider]

I could not get this to work.

Part of the problem may be that with the trigonometric terms, there are likely an infinity of possible results. The ‘z1’ values always plotted on the ‘eqn1’ line, however not at the intersection.

I created the missing ‘‘f1’ function as:

f1 = @(b) eqn1(b(1),b(2)) - eqn2(b(1),b(2));

I moved your Comment here, and deleted my Answer.

.

[Zeyad]

Thanks for taking the time to attempt solving the issue. Do you have any recommendations of what I should do/attempt next? Any ideas are welcome :)

[Star Strider]

My pleasure!

Since I am not certain what the problem is with the functions with respect to finding the intersection, consider using the Global Optimization Toolbox functions, such as ga and others, to see if that approach could solve it. Another possibility in that respect is the fminsearch function. They all use derivative-free algorithms that could avoid whatever the problems are with the gradient-descent algorithms in finding the intersection.

Answer 1

[Alan Weiss]

I don't know if this is helpful, but I suspected that fsolve was having trouble because of the scaling of your variables; everything is so small that fsolve might not be able to tell what is near zero. So I altered the function values by multiplying by 1e3:

f1 = @(er1) 1e3* [u1^2 + w1^2 - ((er1(2)*c1_rb)^2 + (er1(1)*c1_jr)^2 - ...

This gave me the solution

z1 =

  1.000000000000000 + 0.000000000000000i -1.029681858998663 + 0.000000003352866i

Hmm, a little bit complex. What of the function value at the real part of the solution?

pp = real(z1);
f1(pp)

ans =

   1.0e-11 *

  0.014140143187397 - 0.000655983081089i
 -0.510702591327572 + 0.000000132969878i

So the real part of the solution gives a miniscule function value, even after multiplication by 1e3. Is this helpful for you?

Alan Weiss

MATLAB mathematical toolbox documentation

[Zeyad]

First of all, thanks for taking the time to answer. I tried your suggestion of multiplying by 1e3 but I am not getting complex numbers, I am actually getting the same values as the unaltered function, which also changes with the initial guess.