I have a problem with the convergence of fsolve ?

6 views (last 30 days)
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')

Accepted Answer

Matt J
Matt J on 17 Nov 2023
Edited: Matt J on 17 Nov 2023
Your equations look like polynomials of rather large degree, at least 8. The solutions of high order polynomials are known to be unstable.
  3 Comments
Matt J
Matt J on 17 Nov 2023
You could compute the Jacobian at your initial point and then check cond(J).

Sign in to comment.

More Answers (3)

John D'Errico
John D'Errico on 17 Nov 2023
  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.

Torsten
Torsten on 17 Nov 2023
Edited: Torsten on 17 Nov 2023
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.

Matt J
Matt J on 17 Nov 2023
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!