Symbolical differentiation of parametric function

Hello! I am trying to compute the symbolic partial derivative of a parametric function I have defined.

I will report here both main and function

C = [-0.016 1.82e-2 5.5e-3 -2.06e-5 5.91e-7];

G = 6.67259e-20;

R_main = 780e-3;

l = 2;

mass_sys = 5.28e11; %[kg] mass of the system

mass_ratio = 0.0093; % mass ratio

mu = 1 / ((1/mass_ratio) + 1);

m_main = (1 - mu) * mass_sys;

syms x y z

U_main_x = diff(@(x, y, z) U_sph_harm(x, y, z, C, G, m_main, R_main, l), x);

function U = U_sph_harm(x, y, z, C, G, M, R, l)

[lambda, phi, r] = cart2sph(x, y, z);

sum_U = 0;

C_count = 1;

for l_count = 1:l

Plm = legendre(2*l, sin(phi));

for m_count = 0:l_count

sum_U = sum_U + (R/r)^(2*l) * C(C_count) * cos(2 * m * lambda) * Plm(2 * m_count);

C_count = C_count + 1;

end

end

U = G * M / r * (1 + sum_U);

I get the following error

Conversion to logical from sym is not possible.

Error in legendre (line 101)

if ~isreal(x) || ischar(x)|| max(abs(x(:))) > 1

Error in U_sph_harm (line 8)

Plm = legendre(2*l, sin(phi));

Error in main_asteroid (line 92)

U_main_x = diff(@(x, y, z) U_sph_harm(x, y, z, C, G, m_main, R_main, l), x);

Error in sym>funchandle2ref (line 1601)

S = x(S{:});

Error in sym>tomupad (line 1514)

x = funchandle2ref(x);

Error in sym (line 332)

S.s = tomupad(x);

Error in sym/privResolveArgs (line 1102)

argout{k} = sym(arg);

Error in sym/diff (line 29)

args = privResolveArgs(S, invars{:});

Error in main_asteroid (line 92)

U_main_x = diff(@(x, y, z) U_sph_harm(x, y, z, C, G, m_main, R_main, l), x);

Could you please help me sort out why?

Thank you!

Answer 1

[Dyuman Joshi]

The function you are using for Legendre polynomials only accepts numerical inputs.

What you are looking for is legendreP, the symbolic math toolbox function.

However, there are other errors in your code, particularly in the addition line of double for loop -

Plm(2 * m_count)

It's not clear to me what you want do with this line of code.

As you can see below, that Plm is a function of x, y, and z. So, in case you want to find the value of the function, you will need to provide three inputs to it.

C = [-0.016 1.82e-2 5.5e-3 -2.06e-5 5.91e-7];
G = 6.67259e-20;
R_main = 780e-3;
l = 2;
mass_sys = 5.28e11; %[kg] mass of the system
mass_ratio = 0.0093; % mass ratio
mu = 1 / ((1/mass_ratio) + 1);
m_main = (1 - mu) * mass_sys;
syms x y z 
U_main_x = diff(U_sph_harm(x, y, z, C, G, m_main, R_main, l), x)

Plm = 

Array indices must be positive integers or logical values.

Error in indexing (line 936)
R_tilde = builtin('subsref',L_tilde,Idx);

Error in solution>U_sph_harm (line 22)
sum_U = sum_U + (R/r)^(2*l) * C(C_count) * cos(2 * M * lambda) * Plm(2 * m_count);

function U = U_sph_harm(x, y, z, C, G, M, R, l)
[lambda, phi, r] = cart2sph(x, y, z);
sum_U = 0;
C_count = 1;
for l_count = 1:l
    
    %Corrected function
    Plm = legendreP(2*l, sin(phi))
    for m_count = 0:l_count
        %                                                  v correction
        sum_U = sum_U + (R/r)^(2*l) * C(C_count) * cos(2 * M * lambda) * Plm(2 * m_count);
        C_count = C_count + 1
    end
end
U = G * M / r * (1 + sum_U);
end

[Astrid Barzaghi]

U has the following mathematical expression (though the second sum goes from m=0 to l, not from m=1)

The goal is to find the symbolic derivative U_x, U_y and U_z

[Dyuman Joshi]

What is the definition of C?

C has to be a function, it can not be an array, as the limits of the outer sum are from 1 to Inf.

[Astrid Barzaghi]

C is limited because l is limited. This is finding the gravity potential of a body based on spherical harmonics. The model reaches degree 4 so l = 2 (symmetries allow to consider on even contributions to the potential)

For this reason C is a limited and known array that contains [C20 C22 C40 C42 C44].

Thanks for the time your are putting into helping me, hope this clarifies what I am trying to do.

[Dyuman Joshi]

@Astrid Barzaghi

What does the notation P(2l,2m) mean?

or is it P(2l,2m,sin(phi)) ?

[Astrid Barzaghi]

P is the associated legendre function of degree 2l, 2m evaluated in sin(phi).

It turned out that the easiest way to solve the problem was to write myself a function containing the asscoiated legendre functions needed.

The final script is the following

C = [-0.016 1.82e-2 5.5e-3 -2.06e-5 5.91e-7];

G = 6.67259e-20;

R_main = 780e-3;

l = 2;

mass_sys = 5.28e11; %[kg] mass of the system

mass_ratio = 0.0093; % mass ratio

mu = 1 / ((1/mass_ratio) + 1);

m_main = (1 - mu) * mass_sys;

syms x y z

U_main_x = diff(U_sph_harm(x, y, z, C, G, m_main, R_main, l), x);

function U = U_sph_harm(x, y, z, C, G, M, R, l)

[lambda, phi, r] = cart2sph(x, y, z);

sum_U = 0;

C_count = 1;

for l_count = 1:l

for m_count = 0:l_count

l_count_st = num2str(2 * l_count);

m_count_st = num2str(2 * m_count);

lm_st = [l_count_st, m_count_st];

lm = str2double(lm_st);

Plm = ass_legendre(lm);

Plm = Plm(sin(phi));

sum_U = sum_U + (R/r)^(2*l) * C(C_count) * cos(2 * m_count * lambda) * Plm;

C_count = C_count + 1;

end

end

U = G * M / r * (1 + sum_U);

function Plm = ass_legendre(lm)

switch lm

case 20

Plm = @(x) (1/2) * (3 * x^2 - 1);

case 22

Plm = @(x) 3 * (1 - x^2);

case 40

Plm = @(x) (1/8) * (35 * x^4 - 30 * x^2 + 3);

case 42

Plm = @(x) (15/2) * (7 * x^2 - 1) * (1 - x^2);

case 44

Plm = @(x) 105 * (1 - x^2);

end

Thanks again for the help!

[Dyuman Joshi]

Ah, I see, it is the Associate Legendre Polynomial. I was not aware of these, that's why I was confused with the notation.

Well you can generate them in this way -

syms x

AP42 = associate(4,2,x)

AP42(x) = 

AP20 = associate(2,0,x)

AP20(x) = 

function AP = associate(l,m,x)
%Definition taken from Wikipedia
AP(x) = (-1)^m*(1-x^2)^(m/2)*diff(legendreP(l,x),x,m);
end