Asked by Veronica
on 23 Feb 2013

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!

*No products are associated with this question.*

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.

Show 7 older comments

Walter Roberson
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?

Walter Roberson
on 24 Feb 2013

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)

Opportunities for recent engineering grads.

## 0 Comments