How to solve 'Equation solved at initial point' with fsolve?

Hi everyone. I have implemented this code in MATLAB, however I obtain that:

fsolve completed because the vector of function values at the initial point

is near zero as measured by the value of the function tolerance, and

the problem appears regular as measured by the gradient.

How can I solve this problem? Here my code

clc; clear all; close all;

global mu_earth J2 R_earth M_dot_pert Omega_dot_pert

mu_earth = 398600.441;
R_earth = 6378;
M_dot_pert = 1.11337451*10^(-3);
Omega_dot_pert = 1.991063853*10^(-7);
J2=0.0010826269;

fun = @root2d;

%valori di tentativo
x0 = [6858,97.331];

options = optimoptions('fsolve','Display','iter','StepTolerance',1e-12);

x = fsolve(fun,x0,options)

Here the root2d function

function F = root2d(x) 

global mu_earth J2 R_earth M_dot_pert Omega_dot_pert

F(1) = Omega_dot_pert +1.5*J2*((R_earth/x(1))^2)*M_dot_pert*cos(x(2));
F(2) = M_dot_pert - sqrt((mu_earth/x(1)^3))*(1+1.5*J2*(((R_earth/x(1))^2)*(1-1.5*(sin(x(2))^2))));


end

Answer 1

[John D'Errico]

As is often necessary on a problem like this, it helps to look carefully at what you have.

x = sym('x',[1 2]);

mu_earth = 398600.441;
R_earth = 6378;
M_dot_pert = 1.11337451*10^(-3);
Omega_dot_pert = 1.991063853*10^(-7);
J2=0.0010826269;

F(1) = Omega_dot_pert +1.5*J2*((R_earth/x(1))^2)*M_dot_pert*cos(x(2))

F = 

F(2) = M_dot_pert - sqrt((mu_earth/x(1)^3))*(1+1.5*J2*(((R_earth/x(1))^2)*(1-1.5*(sin(x(2))^2))))

F = 

So you have some REALLY HUGE numbers in there. And you are raising large numbers (on the order of 7000, to powers. That alone MAY be a significant problem.

As bad is the idea that you are forming the sin and cos of numbers that are on the order of 100.

Is the number 97 in DEGREES????????? If it is in degrees, then you clearly do not understand that sin and cos assume RADIANs. If you want to work in degrees, then you need to use sind and cosd. In fact, I can be almost positive you think that 97 is a number in degrees, since that is perhaps the most common mistake we ever see on this site.

97 degrees is

deg2rad(97.331)

ans = 1.6987

So 1.6987 RADIANS. Again, if you want to work in degrees, you will need to write the above expressions using sind and cosd. Take your pick. I've done there here. Note that the symbolic toolbox converts degrees to radians, then uses sin and cos.

F(1) = Omega_dot_pert +1.5*J2*((R_earth/x(1))^2)*M_dot_pert*cosd(x(2))

F = 

F(2) = M_dot_pert - sqrt((mu_earth/x(1)^3))*(1+1.5*J2*(((R_earth/x(1))^2)*(1-1.5*(sind(x(2))^2))))

F = 

But at the initial point, now we can look at your objective.

x0 = [6858,97.331];
vpa(subs(F,x,x0),5)

ans = 

And we see here that it results in already very small numbers, near the default tolerance for fsolve.

If I compute the gradient, I'd bet that again, we will see small numbers. You will need to learn to scale your problems, using units that are appropriate, yet will not cause numerical problems. Effectively, you need to learn why it is that astronomers work in light years and mega-parsecs, not in feet or meters.