I would like to understand how to solve this equation with the unknown "a". I admit I'm not an expert with the solve function. I followed the online directions but I continue to receive this error, Someone can help me? Thank you!

syms a j mu D R_earth k pi J_2 
 eqn= j*abs(-2*pi*(2*pi*(a^3/mu)^1/2)/D-(3*pi*J_2*R_earth^2*cos(i))/(a^2*(1-e^2)^2))==k*2*pi;

a(w+1)=solve(eqn,a)

Unable to perform assignment because the left and right sides have a different number of elements.

Error in sym/privsubsasgn (line 1096)
                L_tilde2 = builtin('subsasgn',L_tilde,struct('type','()','subs',{varargin}),R_tilde);

Error in sym/subsasgn (line 933)
            C = privsubsasgn(L,R,inds{:});

Error in part2 (line 37)
a(w+1)=solve(eqn,a)

Answer 1

[Walter Roberson]

e is undefined.

In cos(i) is that intentionally cos(1i) -- cos of the imaginary unit?

We do not know what w is for this purpose.

But imagine that w is a scalar. Then the left side of the assignment would designate space for exactly one output. But is the result of solve() exactly one value?

No! the result of solve() is 5 values. It turns out to be the 5 different roots of a quintic.

Note, by the way, that you are assigning a(w+1) where you just solved for symbolic variable a. There does not seem to be a good reason to do that. Assign to a different variable. If you need the value for a later, then

subs(EXPRESSION, a, solutions_for_a)

[Walter Roberson]

You might get a different number of solutions if you added constraints for your variables. By default the assumption will be that they are complex valued.

[Tommaso Scali]

Ok now I have clearer ideas. Thanks for the indication on the subs function. One last thing before trying if it actually works, I would like to get only the real and positive value, how can I set certain conditions?

[Walter Roberson]

syms a positive

would be one step. However, you have a number of symbolic variables involved, any one of which might be complex valued, so it will not be feasible for MATLAB to figure out which of the roots are going to be positive real. If there were constraints on the variables then it might have a chance .

[Tommaso Scali]

Actually it is a physical problem of an orbit and therefore they are all real numbers. Besides, I already know the value of those variables (I calculate them every time before in the cycle). The problem is I didn't know how to do it without replacing the numbers (which would be impossible given that it is a cycle). I have attached the entire code in order to be clearer.

j=39;
k=3;
for e=0:0.01:0.9 %DOn't put e=1 for singularities 
    
a(1)=((D^2*k^2*mu)/(j^2*4*pi^2))^(1/3); %This is the initial value for a 
   
control=1; %to enter in the cycle

w=1; %index for iteration
    
while abs(control)>1e-10
    
T=2*pi*(a^3/mu)^(1/2);

deltal1=-2*pi*T/D;

deltal2=-(3*pi*J_2*R_earth^2*cos(i))/(a^2*(1-e^2)^2);
 syms a j delta1 delta2 j k 
 eqn= j*abs(delta1+delta2)==k*2*pi;

a(w+1)=solve(eqn,a)

control=a(w+1)-a(w);

w=w+1;
    en
    a=a(end); %Extrapolate the value of a for this ecentricity 

[Walter Roberson]

a(1)=((D^2*k^2*mu)/(j^2*4*pi^2))^(1/3); %This is the initial value for a 

Presumably that is numeric -- though we do not have values for D or mu so it is not certain.

 syms a j delta1 delta2 j k 

This line is equivalent to

 a = sym('a');
 j = sym('j');
 delta1 = sym('delta1');
 delta2 = sym('delta2');
 j = sym('j');
 k = sym('k');
 

Notice this overwrites all of a

[Tommaso Scali]

Yeah everything is numeric. I mean the only unknown in the equation is a. It's an iteration where I'd like to save at every cycle the value and then control if the new value is different or not in order to obtain convergence. I realized that the sym command overwrites all of a, what could I do? Below I attached the entire code also with the definition of data in order to not have misunderstandings anymore, Thank you very much.

%% DATA

format long

mu=398600.441; %mu earth [km^3/s^2]
J_2=1082.631e-6;  %coefficient gravity field earth [-]
R_earth=6378.136; %radius earth [km]
D=86164.1004;  %sideral rotation period [s]
i=deg2rad(70); %For the test of the table 

%% Approach 1
j=39;
k=3;
for e=0:0.01:0.9 %DOn't put e=1 for singularities 
    
a(1)=((D^2*k^2*mu)/(j^2*4*pi^2))^(1/3); %This is the initial value for a 
   
control=1; %to enter in the cycle

w=1; %index for iteration
    
while abs(control)>1e-10
    
T=2*pi*(a^3/mu)^(1/2);

deltal1=-2*pi*T/D;

deltal2=-(3*pi*J_2*R_earth^2*cos(i))/(a^2*(1-e^2)^2);

 syms a j delta1 delta2 j k 
 eqn= j*abs(delta1+delta2)==k*2*pi;

a(w+1)=solve(eqn,a)

control=a(w+1)-a(w);

w=w+1;
    end
    
         
    a=a(end); %Extrapolate the value of a for this ecentricity 

[Walter Roberson]

%% DATA
format long
mu=398600.441; %mu earth [km^3/s^2]
J_2=1082.631e-6;  %coefficient gravity field earth [-]
R_earth=6378.136; %radius earth [km]
D=86164.1004;  %sideral rotation period [s]
i=deg2rad(70); %For the test of the table
%% Approach 1
j=39;
k=3;

syms A positive
I = sym(1i);
PI = sym('pi');

e_vals = 0:0.01:0.9; %DOn't put e=1 for singularities
num_e = length(e_vals);
a_for_e = zeros(1, num_e, 'sym');
for eidx = 1 : num_e
    syms a
    e = e_vals(eidx);
    a(1)=((D^2*k^2*mu)/(j^2*4*PI^2))^(1/3); %This is the initial value for a
    
    control=1; %to enter in the cycle
    w=1; %index for iteration
    
    while abs(control)>1e-10
        
        T=2*PI*(A^3/mu)^(1/2);
        delta1=-2*PI*T/D;
        delta2=-(3*PI*J_2*R_earth^2*cos(I))/(A^2*(1-e^2)^2);
        
        eqn= j*abs(delta1+delta2)==k*2*PI;
        solA = vpasolve(eqn, A, a(w));  %multiple solutions so provide initial estimate
        if isempty(solA)
            a(w+1) = sym(nan);
        else
            a(w+1) = solA;
        end
        control = a(w+1) - a(w);
        w=w+1;
    end
    
    
    a_for_e(eidx) = a(end); %Extrapolate the value of a for this ecentricity
end

plot(e_vals, a_for_e)

Note: there are at least two positive solutions to the equation, so at each point I use the previous solution as the initial value, to try to provide continuity. One of the solutions is in the 7500 range, and the other is in the 1190 range.

[Tommaso Scali]

I finally got it, thank you very much!!

Answer 2

[madhan ravi]

syms a j mu D R_earth k pi J_2 
 eqn= j*abs(-2*pi*(2*pi*(a^3/mu)^1/2)/D-(3*pi*J_2*R_earth^2*cos(i))/(a^2*(1-exp(2))^2))==k*2*pi;
a=solve(eqn,a)