How can I simplify a very complicated derivative (when using simplify does not help much)?

I am working with the following expression:

where

,

,

,

, and

.

My goal is to take the derivative of Π with respect to γ (plugging in and ). Given the complexity of this task, I have written code to do this in MATLAB. Here is my code:

% MATLAB CODE TO CALCULATE THE DERIVATIVE OF PI WITH RESPECT TO GAMMA
% SCRIPT PREPARED ON JANUARY 24, 2024 

% Define symbolic variables
syms gamma theta L q alpha s d1 d2 d3 d4 Pi

% Define d1, d2, d3, and d4
d1 = -L + theta^2*(1 - gamma)^2;
d2 = -L + s*gamma^2 - theta^2*(1 - gamma);
d3 = theta^2 + gamma^alpha - theta^2*gamma^2;
d4 = L + theta^2*(1 - gamma);

% Define s
s = ((gamma^2 * d1 * d3 + 2 * gamma^alpha * (1 - L) * d4) + sqrt(gamma^4 * d1^2 * d3^2 - 4 * gamma^(2 * alpha) * (1 - L) * d1 * d4)) / (2 * gamma^(2 + alpha) * (1 - L));

% Define the expression for Pi
Pi = (1 - q*(1 - L)/d1) * (gamma^alpha - gamma^2 * theta^2) + ...
     (q*(1 - L)*(s * gamma^2 - theta^2 * (1 - gamma) * (2 - gamma)) / (d1 * d2)) * ((1 - s) * gamma^alpha - gamma^2 * theta^2) + ...
     q*(1 - L)/d2 * (gamma^alpha - gamma^2 * theta^2) - ...
     q * (theta^2 + gamma^alpha - gamma^2 * theta^2) * log(q*(1 - L)/d2);

% Calculate the derivative of Pi with respect to gamma
dPi_dgamma = diff(Pi, gamma);

% Display the result
disp(dPi_dgamma);

While this works fine, the result is incredibly long and complicated. I wonder if anyone has any advice on how to simplify things further in MATLAB. I find that using the below code does not do anything at all. Any help or advice would be appreciated! Thanks.

% Simplify the derivative
simplified_dPi_dgamma = simplify(dPi_dgamma)

simplified_dPi_dgamma = 

[Dyuman Joshi]

I doubt there is a short and simple form possible for the obtained expression.

You can try using different name-value arguments of simplify, such as setting the number of steps or the time for simplification, but the results might not be as expceted.

Another thing you can try is changing the output format of the expression by using -

sympref('AbbreviateOutput',false);

This will not abbreviate the output obtained (i.e. there will be no sigmas where terms are grouped).

Answer 1

[Torsten]

You can try

nsteps = 100;
simplified_dPi_dgamma = simplify(dPi_dgamma,'Steps',nsteps) 

[Star Strider]

I would use at least:

nsteps = 500;

It will stop when it can either simplify no further or reaches that limit.

Answer 2

[Morgan]

I'm assuming alpha, theta, q, and L are constants with respect to gamma. This code provides a bit simplier of a final equation and I used pretty() instead of disp() since it reads a bit better in my experience.

% DEFINE CONSTANTS
syms alpha theta q L 

% DEFINE FUNCTIONS OF GAMMA
syms s(gamma) d1(gamma) d2(gamma) d3(gamma) d4(gamma)

% DEFINE PI EXPRESSION
Pi = (1 - q*(1 - L)/d1) * (gamma^alpha - gamma^2 * theta^2) + ...
     (q*(1 - L)*(s * gamma^2 - theta^2 * (1 - gamma) * (2 - gamma)) / (d1 * d2)) * ...
     ((1 - s) * gamma^alpha - gamma^2 * theta^2) + ...
     q*(1 - L)/d2 * (gamma^alpha - gamma^2 * theta^2) - ...
     q * (theta^2 + gamma^alpha - gamma^2 * theta^2) * log(q*(1 - L)/d2);

% CALCULATE DERIVATIVE
dPi_dgamma = diff(Pi, gamma);

% DISPLAY PRETTY RESULT
pretty(dPi_dgamma);

% CONVERT TO LATEX
chr = latex(dPi_dgamma);
clipboard('copy',chr);

Since d1(gamma), d2(gamma), d3(gamma), and d4(gamma) are relatively easy derivatives to take with respect to gamma it shows the result before factoring those terms out. This script also copies the latex equation to your clipboard and with that I get

Not the prettiest, but the equation isn't pretty to begin with...

[Morgan]

With similar method of calculating derivatives, the rest of the derivatives are

[Bryce]

Thanks for the response. Just to clarify--and maybe this reflects a deeper misunderstanding on my part--but why are there still a bunch of derivative symbols ∂ in the result? I don't want that. For context, I am looking to take and then set the result and solve for . Does that make sense?

[Morgan]

Ahhhh in that case add to the end of your code

% SOLVE FOR GAMMA*
gamma_star = solve(dPi_dgamma == 0,gamma)

MATLAB will be unable to find an explicit solution unless you can provide α, and possibly some of the other constants, as an input. If it still won't give a solution take a look at Troubleshoot Equation Solutions from solve Function.

[Morgan]

And to answer your question about why there remains some derivatives in my original answer, I was attempting to give the broadest/simplest possible analytical derivative with fewer "sigma" terms. The actual derivative requires you to substitute the expressions for those remaining partial derivatives (see first comment under this thread) into it. This is equivalent to the derivative you found in your original question, which is likely not any prettier of a solution.

[Bryce]

Thanks. Do you think this is the best way to go about the derivative? In other words, assuming alpha, theta, q, and L are constants with respect to gamma, and then using pretty()? Do you have the code for that? I only seen the code for the first part?

[Torsten]

% DEFINE CONSTANTS
syms alpha theta q L 
% DEFINE FUNCTIONS OF GAMMA
syms s(gamma) d1(gamma) d2(gamma) d3(gamma) d4(gamma)

% Define d1, d3, and d4
d1(gamma) = -L + theta^2*(1 - gamma)^2;
d3(gamma) = theta^2 + gamma^alpha - theta^2*gamma^2;
d4(gamma) = L + theta^2*(1 - gamma);
% Define s
s(gamma) = ((gamma^2 * d1(gamma) * d3(gamma) + 2 * gamma^alpha * (1 - L) * d4(gamma)) + sqrt(gamma^4 * d1(gamma)^2 * d3(gamma)^2 - 4 * gamma^(2 * alpha) * (1 - L) * d1(gamma) * d4(gamma))) / (2 * gamma^(2 + alpha) * (1 - L));
% Define d2
d2(gamma) = -L + s(gamma)*gamma^2 - theta^2*(1 - gamma);
% Compute derivatives with respect to gamma
dd1dgamma = diff(d1,gamma)

dd1dgamma(gamma) = 

dd2dgamma = diff(d2,gamma)

dd2dgamma(gamma) = 

dd3dgamma = diff(d3,gamma)

dd3dgamma(gamma) = 

dd4dgamma = diff(d4,gamma)

dd4dgamma(gamma) = 

dsdgamma = diff(s,gamma)

dsdgamma(gamma) = 

[Morgan]

If you're wanting an analytical expression for to paste into a publication or conference poster/presentation I'd definitely break it up into an expression of the form

then define what you mean by

similar to what I did in previous comments. You might have to do some extra formatting on those equations since latex() and solve() tend to have weird formatting rules. If the intention is to express

you're going to have a very, very, very long expression that will not be very coherent to anyone.

If you just want a code that solves for

given numbers for alpha, theta, q, and L you can use

% DEFINE CONSTANTS
% Replace these with any constant or function 
% that is independent of gamma
alpha = 1;  
theta = 1;
q     = 1;
L     = 1;

% DEFINE SYMBOLIC VARIABLE
syms gamma

% DEFINE d1, d3, and d4
d1 = -L + theta^2*(1 - gamma)^2;;
d3 = theta^2 + gamma^alpha - theta^2*gamma^2;
d4 = L + theta^2*(1 - gamma);

% Define s
s = ((gamma^2 * d1 * d3 + 2 * gamma^alpha * (1 - L) * d4) + ...
    sqrt(gamma^4 * d1^2 * d3^2 - 4 * gamma^(2 * alpha) * (1 - L) * d1 * d4)) / ...
    (2 * gamma^(2 + alpha) * (1 - L));

% DEFINE d2
d2 = -L + s*gamma^2 - theta^2*(1 - gamma);

% DEFINE Pi EXPRESSION
Pi = (1 - q*(1 - L)/d1) * (gamma^alpha - gamma^2 * theta^2) + ...
     (q*(1 - L)*(s * gamma^2 - theta^2 * (1 - gamma) * (2 - gamma)) / (d1 * d2)) * ((1 - s) * gamma^alpha - gamma^2 * theta^2) + ...
     q*(1 - L)/d2 * (gamma^alpha - gamma^2 * theta^2) - ...
     q * (theta^2 + gamma^alpha - gamma^2 * theta^2) * log(q*(1 - L)/d2);

% CALCULATE DERIVATIVE OF Pi
dPi_dgamma = diff(Pi, gamma);

% SOLVE FOR GAMMA*
gamma_star = solve(dPi_dgamma == 0,gamma)

gamma_star = 

[Morgan]

The code I used to get expressions in the first two answers is:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CALCULATE Pi' = f(d1',d2',d3',d4',s')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% DEFINE CONSTANTS
syms alpha theta q L 

% DEFINE FUNCTIONS OF GAMMA
syms s(gamma) d1(gamma) d2(gamma) d3(gamma) d4(gamma)

% DEFINE PI EXPRESSION
Pi = (1 - q*(1 - L)/d1) * (gamma^alpha - gamma^2 * theta^2) + ...
     (q*(1 - L)*(s * gamma^2 - theta^2 * (1 - gamma) * (2 - gamma)) / (d1 * d2)) * ...
     ((1 - s) * gamma^alpha - gamma^2 * theta^2) + ...
     q*(1 - L)/d2 * (gamma^alpha - gamma^2 * theta^2) - ...
     q * (theta^2 + gamma^alpha - gamma^2 * theta^2) * log(q*(1 - L)/d2);

% CALCULATE DERIVATIVE
dPi_dgamma = diff(Pi, gamma);

% DISPLAY PRETTY RESULT
pretty(dPi_dgamma);

                                                                                      alpha        2
/ q (L - 1)     \                        /   q (L - 1) \                   q #2 (gamma      + theta  - #13)   q (alpha #7 - #4) (L - 1)   q #11 #6 (L - 1)   q #11 #2 (L - 1)
| --------- + 1 | (alpha #7 - #4) - q log| - --------- | (alpha #7 - #4) + -------------------------------- - ------------------------- - ---------------- + ----------------
\ d1(gamma)     /                        \   d2(gamma) /                               d2(gamma)                      d2(gamma)                  #9                 #8

                        alpha                                                                  2           2                    2
     q #1 (L - 1) (gamma      #5 + #4 + alpha #7 #12)   q #10 (L - 1) (2 gamma s(gamma) + gamma  #5 - theta  (gamma - 1) - theta  (gamma - 2))   q #10 #1 #6 (L - 1)   q #10 #1 #2 (L - 1)
   + ------------------------------------------------ + -------------------------------------------------------------------------------------- - ------------------- - -------------------
                            #3                                                                    #3                                                 #9 d2(gamma)          d1(gamma) #8

where

              2                 2
   #1 == gamma  s(gamma) - theta  (gamma - 1) (gamma - 2)

            d
   #2 == ------ d2(gamma)
         dgamma

   #3 == d1(gamma) d2(gamma)

                      2
   #4 == 2 gamma theta

            d
   #5 == ------ s(gamma)
         dgamma

            d
   #6 == ------ d1(gamma)
         dgamma

              alpha - 1
   #7 == gamma

                  2
   #8 == d2(gamma)

                  2
   #9 == d1(gamma)

               alpha
   #10 == gamma      #12 + #13

               alpha
   #11 == gamma      - #13

   #12 == s(gamma) - 1

               2      2
   #13 == gamma  theta

% OPTIONAL CODE TO PRINT LATEX EQUATION
dPi_latex = latex(dPi_dgamma);
fprintf('dPi_dgamma = \n%s\n',dPi_latex)

dPi_dgamma = 
\left(\frac{q\,\left(L-1\right)}{d_{1}\left(\gamma \right)}+1\right)\,\left(\alpha \,\gamma ^{\alpha -1}-2\,\gamma \,\theta ^2\right)-q\,\ln\left(-\frac{q\,\left(L-1\right)}{d_{2}\left(\gamma \right)}\right)\,\left(\alpha \,\gamma ^{\alpha -1}-2\,\gamma \,\theta ^2\right)+\frac{q\,\left(\gamma ^\alpha +\theta ^2-\gamma ^2\,\theta ^2\right)\,\frac{\partial }{\partial \gamma } d_{2}\left(\gamma \right)}{d_{2}\left(\gamma \right)}-\frac{q\,\left(\alpha \,\gamma ^{\alpha -1}-2\,\gamma \,\theta ^2\right)\,\left(L-1\right)}{d_{2}\left(\gamma \right)}-\frac{q\,\left(\gamma ^\alpha -\gamma ^2\,\theta ^2\right)\,\left(L-1\right)\,\frac{\partial }{\partial \gamma } d_{1}\left(\gamma \right)}{{d_{1}\left(\gamma \right)}^2}+\frac{q\,\left(\gamma ^\alpha -\gamma ^2\,\theta ^2\right)\,\left(L-1\right)\,\frac{\partial }{\partial \gamma } d_{2}\left(\gamma \right)}{{d_{2}\left(\gamma \right)}^2}+\frac{q\,\left(\gamma ^2\,s\left(\gamma \right)-\theta ^2\,\left(\gamma -1\right)\,\left(\gamma -2\right)\right)\,\left(L-1\right)\,\left(\gamma ^\alpha \,\frac{\partial }{\partial \gamma } s\left(\gamma \right)+2\,\gamma \,\theta ^2+\alpha \,\gamma ^{\alpha -1}\,\left(s\left(\gamma \right)-1\right)\right)}{d_{1}\left(\gamma \right)\,d_{2}\left(\gamma \right)}+\frac{q\,\left(\gamma ^\alpha \,\left(s\left(\gamma \right)-1\right)+\gamma ^2\,\theta ^2\right)\,\left(L-1\right)\,\left(2\,\gamma \,s\left(\gamma \right)+\gamma ^2\,\frac{\partial }{\partial \gamma } s\left(\gamma \right)-\theta ^2\,\left(\gamma -1\right)-\theta ^2\,\left(\gamma -2\right)\right)}{d_{1}\left(\gamma \right)\,d_{2}\left(\gamma \right)}-\frac{q\,\left(\gamma ^\alpha \,\left(s\left(\gamma \right)-1\right)+\gamma ^2\,\theta ^2\right)\,\left(\gamma ^2\,s\left(\gamma \right)-\theta ^2\,\left(\gamma -1\right)\,\left(\gamma -2\right)\right)\,\left(L-1\right)\,\frac{\partial }{\partial \gamma } d_{1}\left(\gamma \right)}{{d_{1}\left(\gamma \right)}^2\,d_{2}\left(\gamma \right)}-\frac{q\,\left(\gamma ^\alpha \,\left(s\left(\gamma \right)-1\right)+\gamma ^2\,\theta ^2\right)\,\left(\gamma ^2\,s\left(\gamma \right)-\theta ^2\,\left(\gamma -1\right)\,\left(\gamma -2\right)\right)\,\left(L-1\right)\,\frac{\partial }{\partial \gamma } d_{2}\left(\gamma \right)}{d_{1}\left(\gamma \right)\,{d_{2}\left(\gamma \right)}^2}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CALCULATE d1', d2', d3', d4', s'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% DEFINE d1, d2, d3, d4, AND s 
d1 = -L + theta^2*(1 - gamma)^2;
d2 = -L + s*gamma^2 - theta^2*(1 - gamma);
d3 = theta^2 + gamma^alpha - theta^2*gamma^2;
d4 = L + theta^2*(1 - gamma);
s = ((gamma^2 * d1 * d3 + 2 * gamma^alpha * (1 - L) * d4) + ...
    sqrt(gamma^4 * d1^2 * d3^2 - 4 * gamma^(2 * alpha) * (1 - L) * d1 * d4)) / ...
    (2 * gamma^(2 + alpha) * (1 - L));

% SOLVE AND DISPLAY FOR d1', d2', d3', d4', AND s'
d1p = diff(d1,gamma);
pretty(d1p);

     2
theta  (2 gamma - 2)

d2p = diff(d2,gamma);
pretty(d2p);

                        2    d                   2
2 gamma s(gamma) + gamma  ------ s(gamma) + theta
                          dgamma

d3p = diff(d3,gamma);
pretty(d3p);

           alpha - 1                2
alpha gamma          - 2 gamma theta

d4p = diff(d4,gamma);
pretty(d4p);

       2
- theta

dsp = diff( s,gamma);
pretty(dsp);

  /        3   2   2        4   2                   2                        2                        4      2         2                2 alpha - 1
  | 4 gamma  #7  #8  + gamma  #7  #4 #8 2 + #2 theta  #7 (L - 1) 4 + #2 theta  #3 #1 (L - 1) 4 - gamma  theta  #7 #3 #8  2 - alpha gamma            #7 #1 (L - 1) 8          alpha      2
- | --------------------------------------------------------------------------------------------------------------------------------------------------------------- + 2 gamma      theta  (L - 1)
  |                                                                                    2
  \                                                      2 sqrt(#5 - 4 #2 #7 (L - theta  (gamma - 1)) (L - 1))

                                                                                \                                                2                                              alpha
                          2              2      2                               |         alpha + 2            (alpha + 2) (gamma  #7 #8 - sqrt(#5 - #2 #7 #1 (L - 1) 4) + gamma      #1 (L - 1) 2)
   - 2 gamma #7 #8 - gamma  #7 #4 + gamma  theta  #3 #8 - alpha #6 #1 (L - 1) 2 |/(2 gamma          (L - 1)) - ------------------------------------------------------------------------------------
                                                                                |                                                                   alpha + 3
                                                                                /                                                            2 gamma          (L - 1)

where

                          2
   #1 == (1 - gamma) theta  + L

              2 alpha
   #2 == gamma

   #3 == 2 gamma - 2

                                 2
   #4 == alpha #6 - 2 gamma theta

              4   2   2
   #5 == gamma  #7  #8

              alpha - 1
   #6 == gamma

                  2            2
   #7 == L - theta  (gamma - 1)

              alpha        2        2      2
   #8 == gamma      + theta  - gamma  theta

% OPTIONAL CODE TO PRINT LATEX EQUATION
d1p_latex = latex(d1p);
fprintf('d1p = \n%s\n',d1p_latex)

d1p = 
\theta ^2\,\left(2\,\gamma -2\right)

d2p_latex = latex(d2p);
fprintf('d2p = \n%s\n',d2p_latex)

d2p = 
2\,\gamma \,s\left(\gamma \right)+\gamma ^2\,\frac{\partial }{\partial \gamma } s\left(\gamma \right)+\theta ^2

d3p_latex = latex(d3p);
fprintf('d3p = \n%s\n',d3p_latex)

d3p = 
\alpha \,\gamma ^{\alpha -1}-2\,\gamma \,\theta ^2

d4p_latex = latex(d4p);
fprintf('d4p = \n%s\n',d4p_latex)

d4p = 
-\theta ^2

dsp_latex = latex(dsp);
fprintf('dsp = \n%s\n',dsp_latex)

dsp = 
-\frac{\frac{4\,\gamma ^3\,{\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)}^2\,{\left(\gamma ^\alpha +\theta ^2-\gamma ^2\,\theta ^2\right)}^2+2\,\gamma ^4\,{\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)}^2\,\left(\alpha \,\gamma ^{\alpha -1}-2\,\gamma \,\theta ^2\right)\,\left(\gamma ^\alpha +\theta ^2-\gamma ^2\,\theta ^2\right)+4\,\gamma ^{2\,\alpha }\,\theta ^2\,\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)\,\left(L-1\right)+4\,\gamma ^{2\,\alpha }\,\theta ^2\,\left(2\,\gamma -2\right)\,\left(L-\theta ^2\,\left(\gamma -1\right)\right)\,\left(L-1\right)-2\,\gamma ^4\,\theta ^2\,\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)\,\left(2\,\gamma -2\right)\,{\left(\gamma ^\alpha +\theta ^2-\gamma ^2\,\theta ^2\right)}^2-8\,\alpha \,\gamma ^{2\,\alpha -1}\,\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)\,\left(L-\theta ^2\,\left(\gamma -1\right)\right)\,\left(L-1\right)}{2\,\sqrt{\gamma ^4\,{\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)}^2\,{\left(\gamma ^\alpha +\theta ^2-\gamma ^2\,\theta ^2\right)}^2-4\,\gamma ^{2\,\alpha }\,\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)\,\left(L-\theta ^2\,\left(\gamma -1\right)\right)\,\left(L-1\right)}}+2\,\gamma ^\alpha \,\theta ^2\,\left(L-1\right)-2\,\gamma \,\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)\,\left(\gamma ^\alpha +\theta ^2-\gamma ^2\,\theta ^2\right)-\gamma ^2\,\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)\,\left(\alpha \,\gamma ^{\alpha -1}-2\,\gamma \,\theta ^2\right)+\gamma ^2\,\theta ^2\,\left(2\,\gamma -2\right)\,\left(\gamma ^\alpha +\theta ^2-\gamma ^2\,\theta ^2\right)-2\,\alpha \,\gamma ^{\alpha -1}\,\left(L-\theta ^2\,\left(\gamma -1\right)\right)\,\left(L-1\right)}{2\,\gamma ^{\alpha +2}\,\left(L-1\right)}-\frac{\left(\alpha +2\right)\,\left(\gamma ^2\,\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)\,\left(\gamma ^\alpha +\theta ^2-\gamma ^2\,\theta ^2\right)-\sqrt{\gamma ^4\,{\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)}^2\,{\left(\gamma ^\alpha +\theta ^2-\gamma ^2\,\theta ^2\right)}^2-4\,\gamma ^{2\,\alpha }\,\left(L-\theta ^2\,{\left(\gamma -1\right)}^2\right)\,\left(L-\theta ^2\,\left(\gamma -1\right)\right)\,\left(L-1\right)}+2\,\gamma ^\alpha \,\left(L-\theta ^2\,\left(\gamma -1\right)\right)\,\left(L-1\right)\right)}{2\,\gamma ^{\alpha +3}\,\left(L-1\right)}