How come an explicit solution cannot be found using solve?

The error is:

Warning: 362 equations in 1 variables. 
> In /Applications/MATLAB_R2012a.app/toolbox/symbolic/symbolic/symengine.p>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
found. 
> 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 )

Thanks!

Answer 1

[Walter Roberson]

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.

[Veronica]

That works! But the solution gives R_GEOID as a vector of 5 with no phi variables. Do you know why this is?

More specifically,

>> S.R_GEOID

ans =

                                            6356734.493731561922064225567129
                                           226610.45796207450510223122358194
                                          -218390.19251969435880732826957049
   - 640.81445634977872349563377206852 + 73434.025222104292667388050155954*i
   - 640.81445634977872349563377206852 - 73434.025222104292667388050155954*i

[Walter Roberson]

Specify the variable you want to solve for,

solve(....., R_GEOID)

the result will be an expression in terms of phi.

However, be advised that the first of the two expressions you are solving for is always true and so is not contributing anything to the situation.

[Veronica]

The equation I was supposed to be using could be simplified to this:

U_f = U_g + U_c;

S = solve((G_M./R_Geoid).*(1 - (C2 + C3 + C4)) == U_f)

But not I'm getting the RootsOf(z) issue for the solved output. I saw your response to another person regarding this issue, but is there a way around this issue? Should I solve using another program?

[Walter Roberson]

Yes, the solution for R_GEOID is a quintic, which MuPAD expresses as a RootOf(). Trying a couple of values, it looks like there would be 3 real solutions for the range you are interested in. Which of the solution(s) are you interested in?

[Veronica]

I was expecting just one solution, so I'm not quite sure. How do you try the values? If I plug in for phi, I still have the z variable. Also, is there a way to plot R_GEOID vs. phi (I can't find a way to convert symbolic to double)?

[Walter Roberson]

The z variable would be within a RootOf(), and would become something MATLAB would solve for. I do not know how many of the roots MATLAB would find.

phis = 0:180;
RG = double(subs(S, phi, phis));
plot(phis, RG)

[Veronica]

The error I get from using double is:

Error using reshape
To RESHAPE the number of elements must
not change.

Error in sym/double (line 705)
            X = reshape(X,siz);

Error in Hyper_HW2_P3 (line 44)
RG = double(subs(S, phi, phis));

So no matter how I do it, I can't seem to convert the sym to double.

[Walter Roberson]

If you use

T = subs(S, phi, phis);

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

[Veronica]

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

[Walter Roberson]

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

double(T)

work?

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));
end

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)