How to solve 'Equation solved at initial point' with fsolve?
6 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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
0 Commenti
Risposte (1)
John D'Errico
il 24 Nov 2022
Modificato: John D'Errico
il 24 Nov 2022
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(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))))
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)
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(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))))
But at the initial point, now we can look at your objective.
x0 = [6858,97.331];
vpa(subs(F,x,x0),5)
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.
0 Commenti
Vedere anche
Categorie
Scopri di più su Loops and Conditional Statements in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!