obvious error in differentiation

Hi, I have a function as

Z = (T-Tp)*lambdam*u + ((p*(h^2)*gammaw)/k2) + ((p^2)*(h^2)*l*alpha*r*gammac)/(D*u*C) + (ro*p*(l^2)*u*gammae)/D;

obviously, this function is composed of four terms. I want to divide Z by h (which is a positive, non-zero variable) and then take the derivative of Z with respect to h and put it equal to 0 and get h. so I write the below code:

dh = diff(Z/h,h);

h1 = solve (dh==0,h);

Now I get

 
h1 = (u*(C*k2*p*(C*D*gammaw*u + alpha*gammac*k2*l*p*r)*(gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D))^(1/2))/(alpha*gammac*k2*l*r*p^2 + C*gammaw*u*D*p)

but if I divide Z by h manually and enter that from the beginning, and then take the derivative, I'll get

h1 = u*((C*k2*(gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D))/(p*(C*D*gammaw*u + alpha*gammac*k2*l*p*r)))^(1/2)

which is consistent with what I get when I do everything manually. Why this happens? what am I missing?

Answer 1

[Torsten]

Did you try

pretty(h1)

and/or

isequal(h1,h2)

( I named the second term h2 ) ?

I think they are equal - at least in h1 you can divide

(C*D*gammaw*u + alpha*gammac*k2*l*p*r)

by

(alpha*gammac*k2*l*r*p^2 + C*gammaw*u*D*p)

to get rid of the quotient .

[Sherwin]

I ran the isequal function and I got 0. If you look at them they are not equal. the first one has a sqrt just on the numerator but in the second one the whole fraction is under sqrt and I am sure that the second one is correct.

[Torsten]

h1 = (u*(C*k2*p*(C*D*gammaw*u + alpha*gammac*k2*l*p*r)*
     (gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D))^(1/2))/
     (alpha*gammac*k2*l*r*p^2 + C*gammaw*u*D*p)
   = u*(C*k2*p)^(1/2)*(C*D*gammaw*u + alpha*gammac*k2*l*p*r)^(1/2)*
     (gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D)^(1/2)/
     ((alpha*gammac*k2*l*r*p + C*gammaw*u*D)*p)
   = u*(C*k2)^(1/2)*p^(1/2)*(C*D*gammaw*u + alpha*gammac*k2*l*p*r)^(1/2)*
     (gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D)^(1/2)/
     ((alpha*gammac*k2*l*r*p + C*gammaw*u*D)*p)
   = u*(C*k2)^(1/2)*(gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D)^(1/2)/
     (p*(alpha*gammac*k2*l*r*p + C*gammaw*u*D))^(1/2)
   = u*((C*k2*(gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D))/
     (p*(alpha*gammac*k2*l*r*p + C*gammaw*u*D)))^(1/2)

[Sherwin]

Thank you so much for your help. However, when I run isequal on your first term:

(u*(C*k2*p*(C*D*gammaw*u + alpha*gammac*k2*l*p*r)*gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D))^(1/2))/(alpha*gammac*k2*l*r*p^2 + C*gammaw*u*D*p)

and your last term:

u*((C*k2*(gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D))/(p*(alpha*gammac*k2*l*r*p + C*gammaw*u*D)))^(1/2)

I still get 0!!! I don't know what I am missing here.

[Steven Lord]

What happens if you use isAlways instead?

syms t
expression1 = sin(t).^2+cos(t).^2;
expression2 = 1;
isAlways(expression1 == expression2)

ans = logical
   1

[Torsten]

@Sherwin

Steven Lord is right: isAlways instead of isequal has to be used.

Take a look at

Test Symbolic Expressions for Equality

under

https://de.mathworks.com/help/symbolic/isequal.html

for the difference between the two.

[Sherwin]

isAlways returns 0 as well. Do you get 1 when you run this?

h1 = (u*(C*k2*p*(C*D*gammaw*u + alpha*gammac*k2*l*p*r)*(gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D))^(1/2))/(alpha*gammac*k2*l*r*p^2 + C*gammaw*u*D*p);

h2 = u*((C*k2*(gammae*p*ro*l^2 + T*lambdam*D - Tp*lambdam*D))/(p*(C*D*gammaw*u + alpha*gammac*k2*l*p*r)))^(1/2);

isAlways(h1==h2)

[Torsten]

Do you get 1 when you run this?

No. Maybe cancellation of

(alpha*gammac*k2*l*r*p^2 + C*gammaw*u*D*p)^(1/2)

makes that the two terms are not always the same.

But I wouldn't agonize over this. Sometimes humans are still more clever than Symbolic Toolboxes.

[Sherwin]

But you see that they are equal, right?

[Torsten]

Under certain assumptions.

When I "show" that h1 = h2, I use e.g. (a*b)^(1/2) = a^(1/2)*b^(1/2). This only holds if a,b are real (not complex) numbers and non-negative.

I suspect that if such assumptions are necessary, MATLAB does not classify two expressions as equal.

But that's a question for the programmers of the symbolic toolbox - I cannot answer it.