I have a problem with the convergence of fsolve ?

I am encountering convergence issues with the Newton-Raphson method implemented using the fsolve function in MATLAB. Despite providing reasonable initial guesses and adjusting solver options, the method fails to converge for my system of equations. I am seeking guidance or suggestions on how to improve the convergence of the Newton-Raphson method in my specific case.

clear all ;
% Define the symbolic variables
syms R1 R2 Omega U
xi = 0.08;
delta = 0.01;
chi = 0.04;
M = 0.01;
a1 = 2.3;
a2 = -18;
lambda = 0.19;
omega1 = 1;
omega2 = 1;
kappa =0.15;
alpha = 0.1;
beta = 0.1;
% Define the equations
EQ1 = (1./16).*beta.^2.*Omega.^2.*R2.^6+((1./2).*lambda.*Omega.^2.*beta-(1./2).*alpha.*Omega.^2.*beta).*R2.^4+(Omega.^4+Omega.^2.*alpha.^2-2.*Omega.^2.*alpha.*lambda+Omega.^2.*lambda.^2-2.*Omega.^2.*omega2+omega2.^2).*R2.^2-kappa.^2.*Omega.^2.*R1.^2;
EQ2 = (9.*M.^2.*Omega.^6.*a2.^2+9.*U.^2.*delta.^2).*R1.^6+(24.*M.^2.*Omega.^4.*U.^2.*a1.*a2-24.*M.*Omega.^4.*U.*xi.*a2-24.*Omega.^2.*U.^2.*delta+24.*U.^2.*delta.*omega1).*R1.^4+(16.*M.^2.*Omega.^2.*U.^4.*a1.^2-32.*M.*Omega.^2.*U.^3.*xi.*a1+16.*Omega.^4.*U.^2+16.*Omega.^2.*U.^2.*xi.^2-32.*Omega.^2.*U.^2.*omega1+16.*U.^2.*omega1.^2).*R1.^2-16.*chi.^2.*U.^2.*R2.^2;
EQ3 = (9./16).*kappa.^2.*Omega.^2.*delta.^2.*R1.^8+(3./2).*kappa.^2.*Omega.^2.*(-Omega.^2+omega1).*delta.*R1.^6+kappa.^2.*Omega.^2.*(-Omega.^2+omega1).^2.*R1.^4-kappa.^2.*Omega.^2.*chi.^2.*R1.^2.*R2.^2+(Omega.^2-omega2).^2.*chi.^2.*R2.^4 ;
% Create a function for the equations
eqns = [EQ1; EQ2; EQ3];

% Define the range of U values
U_values =0:0.001:10;

% Initialize arrays to store solutions
solutions_R1 = [];
solutions_R2 = [];
solutions_Omega = [];

% Create a function handle for fsolve
equations = matlabFunction(eqns, 'Vars', {U, R1, R2, Omega});

% Solve the equations for eachU value
for i = 1:length(U_values)
    U_val =U_values(i);
    
    % Solve the system of equations for the current U value
    initial_guess = [10,10,1.2]; % Initial guess for R1, R2, Omega
    solution = fsolve(@(x) equations(U_val, x(1), x(2), x(3)), initial_guess);
    
    % Store the solutions in arrays
    solutions_R1(i) = solution(1);
    solutions_R2(i) = solution(2);
    solutions_Omega(i) = solution(3);
end

% Plot the solutions
figure(110)
plot(U_values, solutions_R1,'b')
xlabel('\U')
ylabel('R1')
title('R1 vs Beta')

figure(1)
plot(U_values, solutions_R2,'b')
xlabel('\U')
ylabel('R2')
title('R2 vs Beta')

figure(2)
plot(U_values, solutions_Omega,'b')
xlabel('\U')
ylabel('\Omega')
title('Omega vs Beta')

Answer 1

[Matt J]

Your equations look like polynomials of rather large degree, at least 8. The solutions of high order polynomials are known to be unstable.

https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial

[YOUSSEF El MOUSSATI]

Are there other ways to solve the system or study instability?

[Matt J]

You could compute the Jacobian at your initial point and then check cond(J).

[YOUSSEF El MOUSSATI]

Ok thanks for your

Answer 2

[John D'Errico]

1. Solve simpler problems.
2. Choose even better starting values. Yes, you think yours are reasonable, but they are not as good as you want.
3. Solve simpler problems.
4. Solve simpler problems.

I may have repeated myself, but what I have said three times must be true. Ok, everything I said is true.

Your problem will almost certainly have multiple solutions. Remember that a polynomial of degree n will have n solutions, some of them complex, so fsolve will not see them. A system of polynomials will have even more solutions. But also, high degree polynomials will almost certainly pose numerical problems, even working in double precision.

You NEED good starting values. If you are seeing convergence problems, then your starting values are not sufficiently good for this problem.

And due to the issues with "only" approximately 16 digits of dynamic range in a double, expect problems. Ergo my suggestion to find simpler problems to solve.

[YOUSSEF El MOUSSATI]

Thanks

Answer 3

[Torsten]

You could try

solution = [10,10,1.2]; % Initial guess for R1, R2, Omega
% Solve the equations for eachU value
for i = 1:length(U_values)
    U_val = U_values(i);
    
    % Solve the system of equations for the current U value
    initial_guess = solution; % Initial guess for R1, R2, Omega
    solution = fsolve(@(x) equations(U_val, x(1), x(2), x(3)), initial_guess,optimset('TolFun',1e-10,'TolX',1e-10));
    error(i) = norm(equations(U_val, solution(1), solution(2), solution(3)));
    % Store the solutions in arrays
    solutions_R1(i) = solution(1);
    solutions_R2(i) = solution(2);
    solutions_Omega(i) = solution(3);
end

and plot the error against U_values to see if convergence has improved.

But usually for a polynomial system, a specific solution (if any exists) is hard to fix.

[YOUSSEF El MOUSSATI]

OK thanks for your help

Answer 4

[Matt J]

Because you only have 3 unknowns, you could you could do a grid search for the solution(s). Or use contour3 to find the zero level contour.

[YOUSSEF El MOUSSATI]

OK thanks , how I can use contour3 to find the zero level contour?

[Matt J]

As described at the link I gave you.

[YOUSSEF El MOUSSATI]

Ok thank you very much