清空筛选条件
清空筛选条件

I have a problem with the convergence of fsolve ?

3 次查看（过去 30 天）
YOUSSEF El MOUSSATI
YOUSSEF El MOUSSATI 2023-11-17
评论： YOUSSEF El MOUSSATI 2023-11-18
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')

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2023-11-17
编辑：Matt J 2023-11-17
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

更多回答（3 个）

John D'Errico
John D'Errico 2023-11-17
  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 2023-11-17
编辑：Torsten 2023-11-17
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 2023-11-17
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.

类别

在 Help Center 和 File Exchange 中查找有关 Newton-Raphson Method 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by