Advanced Search

Question: Signs of first and second derivatives of bivariate function: fdiff() vs D[]()

Question:Signs of first and second derivatives of bivariate function: fdiff() vs D[]()

MaPal93 170 Maple 2024
differentiation numeric rootof fdiff duplicate-question
May 06 2024
1

I have a complicated bivariate function f(Gamma,rho) that is a RootOf of a quartic. I know that it is strictly positive (one of the four roots at least) for Gamma=0..10 and rho in (-1,+1), with bounds excluded.

I need to find the signs of its first and second derivatives (wrt to Gamma and wrt to rho: 4 derivatives in total).

I encounter numerical issues when I plot3d the derivatives using D[]() vs. fdiff() (numerical function evaluations of the RootOf). I was hoping for the two commands to produce the same output, but they don't it seems. What's going on?

Script:

> 

restart;
_quartic := RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2);
convert(_quartic,radical):
f(Gamma,rho) := simplify(%):

RootOf((8*rho^3+24*rho^2+24*rho+8)*_Z^4+(-12*Gamma*rho^3-12*Gamma*rho^2+12*Gamma*rho+12*Gamma)*_Z^3+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*_Z^2+(4*Gamma*rho^2-4*Gamma)*_Z-Gamma^2*rho^2+2*rho*Gamma^2-Gamma^2)

 (1)


Synthetic representation of derivatives

> 

der1_Gamma := diff(_quartic, Gamma):
der1_rho := diff(_quartic, rho):

Diff('f(Gamma,rho)', Gamma) = collect~(normal(eval(der1_Gamma, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');
Diff('f(Gamma,rho)', rho) = collect~(normal(eval(der1_rho, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');

der2_Gamma := diff(der1_Gamma, Gamma):
der2_rho := diff(der1_rho, rho):

Diff('f(Gamma,rho)', Gamma$2) = collect~(normal(eval(der2_Gamma, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');
Diff('f(Gamma,rho)', rho$2) = collect~(normal(eval(der2_rho, _quartic = 'f(Gamma,rho)')), 'f(Gamma,rho)');

Diff(f(Gamma, rho), Gamma) = -((-6*rho^3-6*rho^2+6*rho+6)*f(Gamma, rho)^3+(5*Gamma*rho^3-5*Gamma*rho^2-5*Gamma*rho+5*Gamma)*f(Gamma, rho)^2+(2*rho^2-2)*f(Gamma, rho)-Gamma*rho^2+2*Gamma*rho-Gamma)/((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)

  

Diff(f(Gamma, rho), rho) = -(1/2)*((24*rho^2+48*rho+24)*f(Gamma, rho)^4+(-36*Gamma*rho^2-24*Gamma*rho+12*Gamma)*f(Gamma, rho)^3+(15*Gamma^2*rho^2-10*Gamma^2*rho-5*Gamma^2-8*rho-8)*f(Gamma, rho)^2+8*Gamma*rho*f(Gamma, rho)-2*rho*Gamma^2+2*Gamma^2)/((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)

  

Diff(f(Gamma, rho), Gamma, Gamma) = ((448*rho^8+1792*rho^7+1792*rho^6-1792*rho^5-4480*rho^4-1792*rho^3+1792*rho^2+1792*rho+448)*f(Gamma, rho)^8+(-1632*Gamma*rho^8-3264*Gamma*rho^7+3264*Gamma*rho^6+9792*Gamma*rho^5-9792*Gamma*rho^3-3264*Gamma*rho^2+3264*Gamma*rho+1632*Gamma)*f(Gamma, rho)^7+(2120*Gamma^2*rho^8-8480*Gamma^2*rho^6-208*rho^7+12720*Gamma^2*rho^4-624*rho^6-208*rho^5-8480*Gamma^2*rho^2+1040*rho^4+1040*rho^3+2120*Gamma^2-208*rho^2-624*rho-208)*f(Gamma, rho)^6+(-1200*Gamma^3*rho^8+2400*Gamma^3*rho^7+2400*Gamma^3*rho^6-7200*Gamma^3*rho^5+640*Gamma*rho^7+640*Gamma*rho^6+7200*Gamma^3*rho^3-1920*Gamma*rho^5-2400*Gamma^3*rho^2-1920*Gamma*rho^4-2400*Gamma^3*rho+1920*Gamma*rho^3+1200*Gamma^3+1920*Gamma*rho^2-640*Gamma*rho-640*Gamma)*f(Gamma, rho)^5+(250*Gamma^4*rho^8-1000*Gamma^4*rho^7+1000*Gamma^4*rho^6+1000*Gamma^4*rho^5-632*Gamma^2*rho^7-2500*Gamma^4*rho^4+632*Gamma^2*rho^6+1000*Gamma^4*rho^3+1896*Gamma^2*rho^5+1000*Gamma^4*rho^2-1896*Gamma^2*rho^4+16*rho^6-1000*Gamma^4*rho-1896*Gamma^2*rho^3+32*rho^5+250*Gamma^4+1896*Gamma^2*rho^2-16*rho^4+632*Gamma^2*rho-64*rho^3-632*Gamma^2-16*rho^2+32*rho+16)*f(Gamma, rho)^4+(240*Gamma^3*rho^7-720*Gamma^3*rho^6+240*Gamma^3*rho^5+1200*Gamma^3*rho^4-32*Gamma*rho^6-1200*Gamma^3*rho^3-240*Gamma^3*rho^2+96*Gamma*rho^4+720*Gamma^3*rho-240*Gamma^3-96*Gamma*rho^2+32*Gamma)*f(Gamma, rho)^3+(-25*Gamma^4*rho^7+125*Gamma^4*rho^6-225*Gamma^4*rho^5+125*Gamma^4*rho^4+125*Gamma^4*rho^3-225*Gamma^4*rho^2+125*Gamma^4*rho-25*Gamma^4)*f(Gamma, rho)^2+(16*Gamma^3*rho^6-64*Gamma^3*rho^5+80*Gamma^3*rho^4-80*Gamma^3*rho^2+64*Gamma^3*rho-16*Gamma^3)*f(Gamma, rho)-5*Gamma^4*rho^6+30*Gamma^4*rho^5-75*Gamma^4*rho^4+100*Gamma^4*rho^3-75*Gamma^4*rho^2+30*Gamma^4*rho-5*Gamma^4)/(((16*rho^2+32*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^2+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^2-10*Gamma^2*rho+5*Gamma^2-4*rho-4)*f(Gamma, rho)+2*Gamma*rho-2*Gamma)*((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)^2)

  

Diff(f(Gamma, rho), rho, rho) = (1/4)*((21504*rho^6+129024*rho^5+322560*rho^4+430080*rho^3+322560*rho^2+129024*rho+21504)*f(Gamma, rho)^10+(-80640*Gamma*rho^6-347136*Gamma*rho^5-526080*Gamma*rho^4-245760*Gamma*rho^3+157440*Gamma*rho^2+199680*Gamma*rho+56064*Gamma)*f(Gamma, rho)^9+(127680*Gamma^2*rho^6+336512*Gamma^2*rho^5+122944*Gamma^2*rho^4-319744*Gamma^2*rho^3-18176*rho^5-246976*Gamma^2*rho^2-90880*rho^4+40576*Gamma^2*rho-181760*rho^3+53696*Gamma^2-181760*rho^2-90880*rho-18176)*f(Gamma, rho)^8+(-110016*Gamma^3*rho^6-112128*Gamma^3*rho^5+191808*Gamma^3*rho^4+172032*Gamma^3*rho^3+57344*Gamma*rho^5-105792*Gamma^3*rho^2+191488*Gamma*rho^4-59904*Gamma^3*rho+192512*Gamma*rho^3+24000*Gamma^3+2048*Gamma*rho^2-94208*Gamma*rho-37888*Gamma)*f(Gamma, rho)^7+(54960*Gamma^4*rho^6-28480*Gamma^4*rho^5-102480*Gamma^4*rho^4+56960*Gamma^4*rho^3-74176*Gamma^2*rho^5+40080*Gamma^4*rho^2-126144*Gamma^2*rho^4-28480*Gamma^4*rho+41088*Gamma^2*rho^3+7440*Gamma^4+138368*Gamma^2*rho^2+5120*rho^4+19776*Gamma^2*rho+20480*rho^3-25536*Gamma^2+30720*rho^2+20480*rho+5120)*f(Gamma, rho)^6+(-15300*Gamma^5*rho^6+30000*Gamma^5*rho^5-2100*Gamma^5*rho^4-24000*Gamma^5*rho^3+50048*Gamma^3*rho^5+14100*Gamma^5*rho^2+5632*Gamma^3*rho^4-6000*Gamma^5*rho-89856*Gamma^3*rho^3+3300*Gamma^5-1024*Gamma^3*rho^2-13056*Gamma*rho^4+39808*Gamma^3*rho-31232*Gamma*rho^3-4608*Gamma^3-15360*Gamma*rho^2+10752*Gamma*rho+7936*Gamma)*f(Gamma, rho)^5+(1875*Gamma^6*rho^6-6250*Gamma^6*rho^5+7125*Gamma^6*rho^4-3500*Gamma^6*rho^3-18316*Gamma^4*rho^5+2125*Gamma^6*rho^2+25540*Gamma^4*rho^4-2250*Gamma^6*rho+12072*Gamma^4*rho^3+875*Gamma^6-26520*Gamma^4*rho^2+12992*Gamma^2*rho^4+6244*Gamma^4*rho+10240*Gamma^2*rho^3+980*Gamma^4-15488*Gamma^2*rho^2-9728*Gamma^2*rho-512*rho^3+3008*Gamma^2-1536*rho^2-1536*rho-512)*f(Gamma, rho)^4+(3320*Gamma^5*rho^5-9240*Gamma^5*rho^4+7600*Gamma^5*rho^3-560*Gamma^5*rho^2-6144*Gamma^3*rho^4-1320*Gamma^5*rho+4992*Gamma^3*rho^3+200*Gamma^5+6784*Gamma^3*rho^2-4992*Gamma^3*rho+1024*Gamma*rho^3-640*Gamma^3+1536*Gamma*rho^2-512*Gamma)*f(Gamma, rho)^3+(-200*Gamma^6*rho^5+800*Gamma^6*rho^4-1200*Gamma^6*rho^3+800*Gamma^6*rho^2+1248*Gamma^4*rho^4-200*Gamma^6*rho-3136*Gamma^4*rho^3+1920*Gamma^4*rho^2+576*Gamma^4*rho-768*Gamma^2*rho^3-608*Gamma^4+768*Gamma^2*rho)*f(Gamma, rho)^2+(-16*Gamma^5*rho^4+192*Gamma^5*rho^3-480*Gamma^5*rho^2+448*Gamma^5*rho+256*Gamma^3*rho^3-144*Gamma^5-384*Gamma^3*rho^2+128*Gamma^3)*f(Gamma, rho)-20*Gamma^6*rho^4+80*Gamma^6*rho^3-120*Gamma^6*rho^2+80*Gamma^6*rho-32*Gamma^4*rho^3-20*Gamma^6+96*Gamma^4*rho^2-96*Gamma^4*rho+32*Gamma^4)/(((16*rho^2+32*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^2+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^2-10*Gamma^2*rho+5*Gamma^2-4*rho-4)*f(Gamma, rho)+2*Gamma*rho-2*Gamma)*((16*rho^3+48*rho^2+48*rho+16)*f(Gamma, rho)^3+(-18*Gamma*rho^3-18*Gamma*rho^2+18*Gamma*rho+18*Gamma)*f(Gamma, rho)^2+(5*Gamma^2*rho^3-5*Gamma^2*rho^2-5*Gamma^2*rho+5*Gamma^2-4*rho^2-8*rho-4)*f(Gamma, rho)+2*Gamma*rho^2-2*Gamma)^2)

 (2)


Signs of derivatives: fdiff (numerical function evaluations of the RootOf) vs. D[]()

> 

restart;
with(plots):
> 

_quartic := RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2):
> 

plot3d(_quartic, Gamma=0..10, rho=-1..+1, labels=[Gamma,rho,Lambda(Gamma,rho)],axesfont=["helvetica","roman",20],labelfont=["helvetica","roman",30]);
 

 

Define it as a f and test it for Gamma=1 and rho=0.5

> 

f := (Gamma,rho) -> RootOf(-8*(rho + 1)^4*_Z^4 + 12*(rho + 1)^3*Gamma*(rho - 1)*_Z^3 - 5*(rho + 1)^2*(-4/5 + Gamma^2*rho^2 + 2*(-2/5 - Gamma^2)*rho + Gamma^2)*_Z^2 - 4*(rho + 1)*Gamma*(rho^2 - 1)*_Z + Gamma^2*(rho + 1)*(rho - 1)^2):
evalf(f(1.0,0.5));

HFloat(0.5110796212870378)

 (3)

Value at zero:

> 

f(0,0):
allvalues(%):
fl := select(is, [allvalues(f(0,0))], positive)[];evalf(%);

(1/2)*2^(1/2)

  

.7071067810

 (4)

Value at infinity (commented out because too slow)

> 

#limit(f(x,y), {x = infinity, y = 0}):
#fh := select(is, [allvalues(%)], positive)[];evalf(%);

Derivative at zero:

> 

allvalues([D[1](f)(0,0)]):
Dfl := %[1][];

-1/4

 (5)

Derivative at a point, evaluated, vs numerical derivative at a point:

> 

D[1](f)(1,0.5):
evalf(%);
fdiff(f(x,y), x, {x = 1.0, y = 0.5});
fdiff(f, [1], [1.0,0.5]);

D[2](f)(1,0.5):
evalf(%);
fdiff(f(x,y), y, {x = 1.0, y = 0.5});
fdiff(f, [2], [1.0,0.5]);

HFloat(-0.05086932918910799)

  

-0.5086932919e-1

  

-0.5086932919e-1

  

HFloat(-0.05166477232109392)

  

-0.5166477232e-1

  

-0.5166477232e-1

 (6)

Can make a function out of fdiff

> 

fDfG := (Gamma,rho) -> fdiff(f, [1], [Gamma,rho]);
fDfr := (Gamma,rho) -> fdiff(f, [2], [Gamma,rho]);

proc (Gamma, rho) options operator, arrow; fdiff(f, [1], [Gamma, rho]) end proc

  

proc (Gamma, rho) options operator, arrow; fdiff(f, [2], [Gamma, rho]) end proc

 (7)

Check for numerical values close to thresholds:

> 

Digits := 15:
evalf('D[1]'(f)(0.1e-8,0.5));fdiff(f, [1], [0.1e-8,0.5]);
evalf('D[1]'(f)(0.1e-7,0.5));fdiff(f, [1], [0.1e-7,0.5]);
evalf('D[1]'(f)(0.1e-5,0.5));fdiff(f, [1], [0.1e-5,0.5]);
evalf('D[1]'(f)(0.00001,0.5));fdiff(f, [1], [0.00001,0.5]);
evalf('D[1]'(f)(0.001,0.5));fdiff(f, [1], [0.001,0.5]);


evalf('D[2]'(f)(1,-0.99));fdiff(f, [2], [1,-0.99]);
evalf('D[2]'(f)(1,-0.97));fdiff(f, [2], [1,-0.97]);
evalf('D[2]'(f)(1,-0.1));fdiff(f, [2], [1,-0.1]);
evalf('D[2]'(f)(1,0.98));fdiff(f, [2], [1,0.98]);
evalf('D[2]'(f)(1,-0.99));fdiff(f, [2], [1,-0.99]);

57735026.8022959

  

57735026.8022959

  

-0.833333329724894e-1

  

-0.833333329724894e-1

  

-0.833332972489415e-1

  

-0.833332972489415e-1

  

-0.833329724894151e-1

  

-0.833329724894151e-1

  

-0.832972489466445e-1

  

-0.832972489466445e-1

  

-223.615892086941

  

-223.615892086941

  

-43.0236130145893

  

-43.0236130145893

  

-.212392503268663

  

-.212392503268663

  

-0.127828146340716e-2

  

-0.127828146340716e-2

  

-223.615892086941

  

-223.615892086941

 (8)

Compare with D (vertical range here to prevent effect of large values from fdiff near zero):

> 

d1G := plot3d([D[1](f), fDfG], 0..10, -0.95..+0.95, view=-0.3..0, color = [red, blue]);
d1r := plot3d([D[2](f), fDfr], 0..10, -0.95..+0.95, color = [red, blue]);

 

 

 

Second derivatives:

> 

evalf('D[1,1]'(f)(1.0,0.5));
fdiff(f, [1, 1], [1.0,0.5]);

evalf('D[2,2]'(f)(1.0,0.5));
fdiff(f, [2, 2], [1.0,0.5]);

fD2fG := (Gamma,rho) -> fdiff(f, [1, 1], [Gamma]);
fD2fr := (Gamma,rho) -> fdiff(f, [2, 2], [Gamma]);

0.266607527050519e-1

  

0.266607527050519e-1

  

.151600577769391

  

.151600577769391

  

proc (Gamma, rho) options operator, arrow; fdiff(f, [1, 1], [Gamma]) end proc

  

proc (Gamma, rho) options operator, arrow; fdiff(f, [2, 2], [Gamma]) end proc

 (9)
> 

d2G:= plot3d([D[1,1](f), fD2fG], 0..10, -0.9..+0.9, color = [red, blue]);
d2r:= plot3d([D[2,2](f), fD2fr], 0..10, -0.9..+0.9, color = [red, blue]);
 

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

  

 

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

  

 
> 

d1d2G := plot3d([fDfG, fD2fG], 0.1e-6 .. 10, -0.98 .. +0.98, axesfont=["helvetica","roman",20],labelfont=["helvetica","roman",30], size=[1000,1000]);
d1d2r := plot3d([fDfr, fD2fr], 0.1e-6 .. 10, -0.98 .. +0.98, axesfont=["helvetica","roman",20],labelfont=["helvetica","roman",30], size=[1000,1000]);

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

  

 

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

  

 
 

NULL

Download signs_derivatves_bivariate.mw
Please Wait...