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

New to MATLAB?

How come an explicit solution cannot be found using solve?

Asked by Veronica

Veronica (view profile)

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 )




Veronica (view profile)


No products are associated with this question.

1 Answer

Answer by Walter Roberson

Walter Roberson (view profile)

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.


Walter Roberson

Walter Roberson (view profile)

on 24 Feb 2013

If you use

T = subs(S, phi, phis);

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


Veronica (view profile)

on 24 Feb 2013

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

Walter Roberson

Walter Roberson (view profile)

on 24 Feb 2013

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)

Walter Roberson

Walter Roberson (view profile)

Contact us