MATLAB Answers


How come an explicit solution cannot be found using solve?

Asked by Veronica
on 23 Feb 2013

The error is:

Warning: 362 equations in 1 variables. 
> In /Applications/>symengine at 54
  In mupadengine.mupadengine>mupadengine.evalin at 97
  In mupadengine.mupadengine>mupadengine.feval at 150
  In solve at 160
  In Hyper_HW2_P3 at 39 
Warning: Explicit solution could not be
> In solve at 169
  In Hyper_HW2_P3 at 39 

The code is:

% Given parameters of geoid for Earth
U = 6.263685953707e7;  % potential [m^2/s^2]
GM = 3.986005e14;  % gravitational constant [m^3/s^2]
J_2 = 1.08263e-3;  % Jeffrey constants
J_3 = 2.532153e-7;
J_4 = 1.6109876e-7;
w = 7.292115147e-5;  % angular acceleration [rad/s]
R_e = 6.378135e6;  % equatorial radius [m]
R_p = 6.3567506e6;  % polar radius [m]
phi = [0:180];  % angles of colatitude [deg]
% Relations for Legendre Polynomials for scalar potential equation
P_2 = (1/2).*(3.*cosd(phi).^2-1);
P_3 = (1/2).*(5.*cosd(phi).^3 - 3.*cosd(phi));
P_4 = (1/8).*(35.*cosd(phi).^3 - 30.*cosd(phi).^2 + 3);
% Determine R_RE from ellipse equation
R_RE = (R_e*R_p)./(sqrt(R_p^2.*sind(phi).^2 + R_e^2.*cosd(phi).^2));
% Determine potential from the centrifugal force to include the effects of
% rotation for the geoid
syms R_GEOID
U_c = -(1/2).*w^2.*R_GEOID^2.*sind(phi).^2;
% Analytically solve for R_GEOID from potential function
P2 = (R_e/R_GEOID)^2*J_2*P_2;
P3 = (R_e/R_GEOID)^3*J_3*P_3;
P4 = (R_e/R_GEOID)^4*J_4*P_4;
S = solve(-(1/2).*w^2.*R_GEOID.^2.*sind(phi).^2 == U_c,...
    (GM./R_GEOID).*(1 - (P2 + P3 + P4)) == U )




No products are associated with this question.

1 Answer

Answer by Walter Roberson
on 23 Feb 2013
 Accepted answer

When you solve() a matrix, it tries to solve all of the entries simultaneously. It does not consider the elements to be a set of related problems that are to be solved for individually with an array of different results.

What you should do is leave phi as a symbolic parameter instead of making it numeric, and solve() for that, giving you a formula (or set of formulae) that is parameterized on phi. You can then subs() the actual numeric phi values in for the list of solutions.


If you use

T = subs(S, phi, phis);

then what does size(T) give? Also, what is size(phis) for this?

size(T) and size(phis) both give 1 181.

Hmmmm -- and if you do that assignment like that, does



Have a look at a few of the elements of T, and double() them individually to see if double() is trying to return multiple values. If so then

for K = 1 : length(T)
  RG{K} = double(T(K));

If you want to select down to real roots you can then

rRG = cellfun(@(C) C(imag(C)==0), 'Uniform', 0);

And if you are lucky exactly one element per index will remain and then you would

rRGv = cell2mat(rRG);
plot(phis, rRGv)

My analysis so far suggests that 1 of the roots is always real (for real values of b), and that the other 4 roots have non-zero imaginary parts that approach 0 in the limiting case as b approaches infinity (i.e., only one real root until b = infinity)

Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

MATLAB Academy

New to MATLAB?

Learn MATLAB today!