## How come an explicit solution cannot be found using solve?

on 23 Feb 2013

### Walter Roberson (view profile)

The error is:

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

## Products

No products are associated with this question.

### Walter Roberson (view profile)

on 23 Feb 2013

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

### 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

```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)

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