Entropy Gradient Flows

Needs["VectorAnalysis`"]

Needs["VectorFieldPlots`"]

Weibull

Weibull Entropy Using μ, γ

Clear[W]
W[γ_,β_]:=  (β/γ) ((((x/γ)))^(β-1))   (E^(-((((x/γ)))^β)))

FullSimplify[(∫_0^∞)x W[γ,β] \[DifferentialD]x,{γ>0,β>0}]

FullSimplify[(∫_0^∞)(x^2) W[γ,β] \[DifferentialD]x,{γ>0,β>0}] -(((γ Gamma[1+(1/β)]))^2)

μ= γ Gamma[1+(1/β)];
σ=γ (Sqrt[(Gamma[((2+β)/β)]-((Gamma[1+(1/β)])^2))]);

(γ (Sqrt[(Gamma[((2+β)/β)]-((Gamma[1+(1/β)])^2))]))/ (γ Gamma[1+(1/β)])

((Sqrt[-((Gamma[1+(1/β)])^2)+Gamma[((2+β)/β)]])/(Gamma[1+(1/β)]))/.β->0.577

((Sqrt[-((Gamma[1+(1/β)])^2)+Gamma[((2+β)/β)]])/(Gamma[1+(1/β)]))/.β->1/n

Table[((Sqrt[Gamma[(2+(1/n)) n]-((Gamma[1+n])^2)])/(Gamma[1+n])),{n,1,10}]

Table[(Sqrt[(((2 n)!)/(((n!))^2))-1]),{n,1,10}]

Log[(((E^(-((((x/γ)))^β))) β ((((x/γ)))^β))/x)]//PowerExpand

FullSimplify[-(x^β) (γ^(-β))-Log[x]+Log[β]+β (Log[x]-Log[γ])( (β/γ) ((((x/γ)))^(β-1))   (E^(-((((x/γ)))^β))) ),{γ>0,β>0}]

Limit[(∫_z^∞) -(-Log[x]+Log[β]+((((x/γ)))^β) (-1+(((E^(-((((x/γ)))^β))) (β^2) Log[(x/γ)])/x)))\[DifferentialD]x,z->0]//FullSimplify

FullSimplify[f log(f)]

FullSimplify[-(∫_0^∞)(((E^(-(((x θ))^β))) β (((x θ))^β) ((-(((x θ))^β)+log(β)+β log(x θ))-log(x)))/x)\[DifferentialD]x,{(Subscript[α, 1])>0,(Subscript[α, 2])>0}]

-log(β)-log(θ)-(ℽ/β)+ℽ+1/.θ->(1/ξ)

Plot3D[1+EulerGamma-(EulerGamma/β)-Log[β]-Log[(1/γ)],{γ,0,5},{β,0,4}]

WGrad=GradientFieldPlot[1+EulerGamma-(EulerGamma/β)-Log[β]-Log[(1/γ)],{γ,0.01,5},{β,0.001,4},Frame->True,GridLines->Automatic]

WEntCont=ContourPlot[1+EulerGamma-(EulerGamma/β)-Log[β]-Log[(1/γ)],{γ,0.01,5},{β,0.001,4},Contours->50,ColorFunction->(Hue[#]&), ClippingStyle->Automatic]

WcontGrad=Show[WEntCont,WGrad]

Gradient flow curves in ξ,β

S[ξ_,β_]:=1+EulerGamma-(EulerGamma/β)-Log[β]-Log[(1/ξ)]

Plot[1+EulerGamma-(EulerGamma/β)-Log[β]-Log[(1/110)],{β,0,4},GridLines->Automatic,Frame->True]

S[ξ,β]

{(∂_ξ) S[ξ,β],(∂_β) S[ξ,β]}//FullSimplify

{(1/ξ),((EulerGamma-β)/(β^2))}/.{ξ->x[t],β->y[t]}

NDSolve[{x'[t]==(1/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==.1,y[0]==4},{x[t],y[t]},{t,-1,40}]
w0140=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.004,30},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5.1},{0,4.1}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==.1,y[0]==3},{x[t],y[t]},{t,-10,40}]
w0130=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.00,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5.1},{0,4.1}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==.1,y[0]==2},{x[t],y[t]},{t,-10,40}]
w0120=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.00,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5.1},{0,4.1}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==.1,y[0]==1},{x[t],y[t]},{t,-10,40}]
w0110=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.00,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5.1},{0,4.1}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==.1,y[0]==.5},{x[t],y[t]},{t,-10,40}]
w0105=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.00,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5.1},{0,4.1}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==.1,y[0]==.1},{x[t],y[t]},{t,-10,40}]
w0101=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.00,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5.1},{0,4.1}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==1,y[0]==.1},{x[t],y[t]},{t,-10,40}]
w1001=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.00,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5.1},{0,4.1}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==2,y[0]==.1},{x[t],y[t]},{t,-10,40}]
w2001=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.00,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5.1},{0,4.1}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==3,y[0]==.1},{x[t],y[t]},{t,-10,40}]
w3001=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.00,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5.1},{0,4.1}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==4,y[0]==.1},{x[t],y[t]},{t,-10,40}]
w4001=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.00,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5.1},{0,4.1}}]

WintC=Show[w0130,w0140,w0120,w0110,w0105,w0101,w1001,w2001,w3001,w4001]
Show[WcontGrad,WintC]

Weibull entropy using β, θ

Clear[W]

W[θ_,β_]:=  β  θ (((θ x))^(β-1))   (E^(-(((θ x))^β)))

FullSimplify[(∫_0^∞)x W[θ,β] \[DifferentialD]x,{θ>0,β>0}]

FullSimplify[(∫_0^∞)(x^2) W[θ,β] \[DifferentialD]x,{θ>0,β>0}] -(((((Gamma[1+(1/β)])/θ)))^2)

μ= ((Gamma[1+(1/β)])/θ);
σ= (Sqrt[(-(((Gamma[1+(1/β)])^2)/(θ^2))+((Gamma[((2+β)/β)])/(θ^2)))]);

FullSimplify[f log(f)]

FullSimplify[-(∫_0^∞)(((E^(-(((x θ))^β))) β (((x θ))^β) ((-(((x θ))^β)+log(β)+β log(x θ))-log(x)))/x)\[DifferentialD]x,{(Subscript[α, 1])>0,(Subscript[α, 2])>0}]

-log(β)-log(θ)-(ℽ/β)+ℽ+1/.θ->(1/ξ)

SN=-log(β)-log(θ)-(EulerGamma/β)+EulerGamma+1;

FullSimplify[-(∫_0^∞)(((E^(-(((x  θ))^β))) β (((x θ))^β) ((-(((x  θ))^β)+log(β)+β log(x  θ))-log(x)))/x)\[DifferentialD]x,{β>0,θ>0}]

Plot3D[1+EulerGamma-(EulerGamma/β)-Log[β]-Log[θ], {θ,0,5},{β,0,3}]

Unit mean has θ= Gamma[1+(1/β)]

1+EulerGamma-(EulerGamma/β)-Log[β]-Log[θ] /. θ-> Gamma[1+(1/β)]

Plot[1+EulerGamma-(EulerGamma/β)-Log[β]-Log[Gamma[1+(1/β)]],{β,0,30}]

Gradient flow curves in θ, β

S[θ_,β_]:=1+EulerGamma-(EulerGamma/β)-Log[β]-Log[θ]

{(∂_θ) S[θ,β],(∂_β) S[θ,β]}//FullSimplify

Clear[x,y,t,θ,β]

{-(1/θ),((EulerGamma-β)/(β^2))} /. {θ-> x[t] ,  β-> y[t]}

NDSolve[{x'[t]==((-1)/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==2,y[0]==1},{x[t],y[t]},{t,-100,10}]
p2010=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-100,1.99},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015]]

NDSolve[{x'[t]==((-1)/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==.1,y[0]==4},{x[t],y[t]},{t,-150,10}]
w0140=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-100,.0049},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,10},{0,10}}]

NDSolve[{x'[t]==((-1)/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==.1,y[0]==2},{x[t],y[t]},{t,-150,10}]
w0120=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-100,.0049},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,10},{0,10}}]

NDSolve[{x'[t]==((-1)/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==.1,y[0]==1},{x[t],y[t]},{t,-150,10}]
w0110=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-100,.0049},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,10},{0,10}}]

NDSolve[{x'[t]==((-1)/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==.1,y[0]==0.5},{x[t],y[t]},{t,-150,10}]
w0105 =ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,-0.2566172265186751`,0.004999985917545043`},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,10},{0,10}}]

NDSolve[{x'[t]==((-1)/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==2,y[0]==0.1},{x[t],y[t]},{t,-150,10}]
w2001 =ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,2},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,10},{0,10}}]

NDSolve[{x'[t]==((-1)/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==4,y[0]==0.1},{x[t],y[t]},{t,-150,10}]
w4001 =ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,8},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,10},{0,10}}]

NDSolve[{x'[t]==((-1)/(x[t])),y'[t]==((EulerGamma-y[t])/((y[t])^2)),x[0]==6,y[0]==0.1},{x[t],y[t]},{t,-150,10}]
w6001 =ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,10},{0,10}}]

Show[w6001,w4001,w2001,w0110,w0120,w0140]

Gamma

Gamma in terms of mean μ and sd σ

((((κ/μ)))^κ) (((E^(-x κ/μ)) (x^(-1+κ)) )/( Gamma[κ]))/.κ->((((μ/σ)))^2)//FullSimplify

(∫_0^∞)(((E^(-((x μ)/(σ^2)))) (x^(-1+((μ^2)/(σ^2)))) ((((μ/(σ^2))))^((μ^2)/(σ^2))))/(Gamma[((μ^2)/(σ^2))]))\[DifferentialD]x/;σ>0&&μ>0//FullSimplify

(∫_0^∞)(((E^(-((x μ)/(σ^2)))) (x^(-1+((μ^2)/(σ^2)))) ((((μ/(σ^2))))^((μ^2)/(σ^2))))/(Gamma[((μ^2)/(σ^2))]))x \[DifferentialD]x/;σ>0&&μ>0//FullSimplify

(1/(Gamma[((μ^2)/(σ^2))]))((((μ/(σ^2))))^((μ^2)/(σ^2)))μ ((((μ/(σ^2))))^(-((μ^2)/(σ^2)))) Gamma[((μ^2)/(σ^2))]//FullSimplify

(∫_0^∞)(((E^(-((x μ)/(σ^2)))) (x^(-1+((μ^2)/(σ^2)))) ((((μ/(σ^2))))^((μ^2)/(σ^2))))/(Gamma[((μ^2)/(σ^2))])) (((x-μ))^2) \[DifferentialD]x  /;σ>0&&μ>0//FullSimplify

FullSimplify[(1/(Gamma[((μ^2)/(σ^2))]))((((μ/(σ^2))))^((μ^2)/(σ^2))) ((((μ/(σ^2))))^(-2-((μ^2)/(σ^2)))) Gamma[1+((μ^2)/(σ^2))] ]

Entropy

From calculation in μ, κ we convert to μ, σ

κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ]/. κ->((((μ/σ)))^2)

S[μ_,σ_]:=((μ^2)/(σ^2))-Log[(μ/(σ^2))]+Log[Gamma[((μ^2)/(σ^2))]]-(-1+((μ^2)/(σ^2))) PolyGamma[0,((μ^2)/(σ^2))]

S[μ,σ]//FullSimplify

((μ^2)/(σ^2))-Log[(μ/(σ^2))]+Log[Gamma[((μ^2)/(σ^2))]]+(1-((μ^2)/(σ^2))) PolyGamma[0,((μ^2)/(σ^2))]/. σ->(μ/(Sqrt[κ]))//FullSimplify

{(∂_μ) S[μ,σ],(∂_σ) S[μ,σ]}

{(∂_μ) (κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ]),(∂_κ) (κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ])}//FullSimplify

{(1/μ),-(((-1+κ) (-1+κ PolyGamma[1,κ]))/κ)}//TeXForm

Solve[(∂_κ) (κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ])==0,{μ,κ}]

Solve[{(∂_μ)(κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ])==0,(∂_κ)(κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ])==0},{μ,κ}]

Reduce[{(∂_μ) S[μ,σ]==0,(∂_σ) S[μ,σ]==0},{μ,σ}]

Solve[(∂_μ) S[μ,σ]==0,σ]

Solve[(∂_σ) S[μ,σ]==0,σ]

Plot[S[3,σ],{σ,.01,5}, PlotRange->All]

S[μ_,σ_]:=((μ^2)/(σ^2))-Log[(μ/(σ^2))]+Log[Gamma[((μ^2)/(σ^2))]]+(1-((μ^2)/(σ^2))) PolyGamma[0,((μ^2)/(σ^2))]

S[μ,σ]

Plot[{S[5,σ],S[6,σ],S[7,σ],S[8,σ],S[9,σ],S[10,σ]},{σ,4,13},Frame->True,GridLines->Automatic]

GGrad=GradientFieldPlot[((μ^2)/(σ^2))-Log[(μ/(σ^2))]+Log[Gamma[((μ^2)/(σ^2))]]-(-1+((μ^2)/(σ^2))) PolyGamma[0,((μ^2)/(σ^2))], {μ,0.001,10},{σ,0.001,8},Frame->True,GridLines->Automatic]

Limit[{(∂_μ) S[μ,σ],(∂_σ) S[μ,σ]},μ->σ]//FullSimplify

{(∂_μ) S[μ,σ],(∂_σ) S[μ,σ]}//FullSimplify

{-(1/μ)+((2 μ)/(σ^2))+((2 μ (-(μ^2)+(σ^2)) PolyGamma[1,((μ^2)/(σ^2))])/(σ^4)),((2 (-(μ^2)+(σ^2)) ((σ^2)-(μ^2) PolyGamma[1,((μ^2)/(σ^2))]))/(σ^5))}/.{μ->x[t],σ->y[t]}

NDSolve[{x'[t]==-(1/(x[t]))+((2 x[t])/((y[t])^2))+((2 PolyGamma[1,(((x[t])^2)/((y[t])^2))] x[t] (-((x[t])^2)+((y[t])^2)))/((y[t])^4)),y'[t]==((2 (-((x[t])^2)+((y[t])^2)) (-PolyGamma[1,(((x[t])^2)/((y[t])^2))] ((x[t])^2)+((y[t])^2)))/((y[t])^5)),x[0]==.01,y[0]==.01},{x[t],y[t]},{t,0,50}]
p0101=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,50},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015]]

NDSolve[{x'[t]==-(1/(x[t]))+((2 x[t])/((y[t])^2))+((2 PolyGamma[1,(((x[t])^2)/((y[t])^2))] x[t] (-((x[t])^2)+((y[t])^2)))/((y[t])^4)),y'[t]==((2 (-((x[t])^2)+((y[t])^2)) (-PolyGamma[1,(((x[t])^2)/((y[t])^2))] ((x[t])^2)+((y[t])^2)))/((y[t])^5)),x[0]==.01,y[0]==1},{x[t],y[t]},{t,-1,5}];
p1010=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,5},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015]];

NDSolve[{x'[t]==-(1/(x[t]))+((2 x[t])/((y[t])^2))+((2 PolyGamma[1,(((x[t])^2)/((y[t])^2))] x[t] (-((x[t])^2)+((y[t])^2)))/((y[t])^4)),y'[t]==((2 (-((x[t])^2)+((y[t])^2)) (-PolyGamma[1,(((x[t])^2)/((y[t])^2))] ((x[t])^2)+((y[t])^2)))/((y[t])^5)),x[0]==.01,y[0]==5},{x[t],y[t]},{t,0,100}];
p1050=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,100},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015]];

NDSolve[{x'[t]==-(1/(x[t]))+((2 x[t])/((y[t])^2))+((2 PolyGamma[1,(((x[t])^2)/((y[t])^2))] x[t] (-((x[t])^2)+((y[t])^2)))/((y[t])^4)),y'[t]==((2 (-((x[t])^2)+((y[t])^2)) (-PolyGamma[1,(((x[t])^2)/((y[t])^2))] ((x[t])^2)+((y[t])^2)))/((y[t])^5)),x[0]==5,y[0]==.01},{x[t],y[t]},{t,0,200}];
p5010=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,200},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015]];

NDSolve[{x'[t]==-(1/(x[t]))+((2 x[t])/((y[t])^2))+((2 PolyGamma[1,(((x[t])^2)/((y[t])^2))] x[t] (-((x[t])^2)+((y[t])^2)))/((y[t])^4)),y'[t]==((2 (-((x[t])^2)+((y[t])^2)) (-PolyGamma[1,(((x[t])^2)/((y[t])^2))] ((x[t])^2)+((y[t])^2)))/((y[t])^5)),x[0]==7,y[0]==.01},{x[t],y[t]},{t,0,600}];
p7010=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,600},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015]];

NDSolve[{x'[t]==-(1/(x[t]))+((2 x[t])/((y[t])^2))+((2 PolyGamma[1,(((x[t])^2)/((y[t])^2))] x[t] (-((x[t])^2)+((y[t])^2)))/((y[t])^4)),y'[t]==((2 (-((x[t])^2)+((y[t])^2)) (-PolyGamma[1,(((x[t])^2)/((y[t])^2))] ((x[t])^2)+((y[t])^2)))/((y[t])^5)),x[0]==3,y[0]==.01},{x[t],y[t]},{t,0,600}];
p3010=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,600},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015]];

NDSolve[{x'[t]==-(1/(x[t]))+((2 x[t])/((y[t])^2))+((2 PolyGamma[1,(((x[t])^2)/((y[t])^2))] x[t] (-((x[t])^2)+((y[t])^2)))/((y[t])^4)),y'[t]==((2 (-((x[t])^2)+((y[t])^2)) (-PolyGamma[1,(((x[t])^2)/((y[t])^2))] ((x[t])^2)+((y[t])^2)))/((y[t])^5)),x[0]==.01,y[0]==3},{x[t],y[t]},{t,0,600}];
p1030=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,600},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015]];

NDSolve[{x'[t]==-(1/(x[t]))+((2 x[t])/((y[t])^2))+((2 PolyGamma[1,(((x[t])^2)/((y[t])^2))] x[t] (-((x[t])^2)+((y[t])^2)))/((y[t])^4)),y'[t]==((2 (-((x[t])^2)+((y[t])^2)) (-PolyGamma[1,(((x[t])^2)/((y[t])^2))] ((x[t])^2)+((y[t])^2)))/((y[t])^5)),x[0]==.01,y[0]==7},{x[t],y[t]},{t,0,600}];
p1070=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,600},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015]];

NDSolve[{x'[t]==-(1/(x[t]))+((2 x[t])/((y[t])^2))+((2 PolyGamma[1,(((x[t])^2)/((y[t])^2))] x[t] (-((x[t])^2)+((y[t])^2)))/((y[t])^4)),y'[t]==((2 (-((x[t])^2)+((y[t])^2)) (-PolyGamma[1,(((x[t])^2)/((y[t])^2))] ((x[t])^2)+((y[t])^2)))/((y[t])^5)),x[0]==1,y[0]==.01},{x[t],y[t]},{t,0,600}];
p1001=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,600},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015]];

GamEntFlow= Show[p0101,p1010,p1050,p7010,p5010,p3010,p1030,p1070,p1001,PlotRange->{{-0.1,10},{0,8}}]

GamEntCont=ContourPlot[((μ^2)/(σ^2))-Log[(μ/(σ^2))]+Log[Gamma[((μ^2)/(σ^2))]]+(1-((μ^2)/(σ^2))) PolyGamma[0,((μ^2)/(σ^2))],{μ,0.01,10},{σ,0.01,8},Contours->50,ColorFunction->(Hue[#]&), ClippingStyle->Automatic]

Show[GamEntCont,GamEntFlow]

Show[GamEntCont,GGrad,GamEntFlow]

ESurf=Plot3D[((μ^2)/(σ^2))-Log[(μ/(σ^2))]+Log[Gamma[((μ^2)/(σ^2))]]+(1-((μ^2)/(σ^2))) PolyGamma[0,((μ^2)/(σ^2))],{μ,0.01,10},{σ,0.01,8}]

GraphicsArray[{ESurfMuK,ESurf } ]

Gamma entropy in μ, σ

Log[(((E^(-((x μ)/(σ^2)))) (x^(-1+((μ^2)/(σ^2)))) ((((μ/(σ^2))))^((μ^2)/(σ^2))))/(Gamma[((μ^2)/(σ^2))]))](((E^(-((x μ)/(σ^2)))) (x^(-1+((μ^2)/(σ^2)))) ((((μ/(σ^2))))^((μ^2)/(σ^2))))/(Gamma[((μ^2)/(σ^2))]))//PowerExpand//FullSimplify

(∫_0^∞)-(1/(Gamma[((μ^2)/(σ^2))]))(E^(-((x μ)/(σ^2)))) (x^(-1+((μ^2)/(σ^2)))) (μ^((μ^2)/(σ^2))) (σ^(-2-((2 (μ^2))/(σ^2)))) ((-(μ^2)+(σ^2)) Log[x]+μ (x-μ Log[μ]+2 μ Log[σ])+(σ^2) Log[Gamma[((μ^2)/(σ^2))]])\[DifferentialD]x//FullSimplify

-(1/(Gamma[((μ^2)/(σ^2))]))(μ^((μ^2)/(σ^2))) (σ^(-2-((2 (μ^2))/(σ^2))))((((μ/(σ^2))))^(-((μ^2)/(σ^2)))) Gamma[((μ^2)/(σ^2))] ((μ^2)-(μ^2) Log[μ]+((μ^2)-(σ^2)) Log[(μ/(σ^2))]+2 (μ^2) Log[σ]+(σ^2) Log[Gamma[((μ^2)/(σ^2))]]+(-(μ^2)+(σ^2)) PolyGamma[0,((μ^2)/(σ^2))])//FullSimplify

(∂_μ) (-(μ^((μ^2)/(σ^2))) ((((μ/(σ^2))))^(-((μ^2)/(σ^2)))) (σ^(-2-((2 (μ^2))/(σ^2)))) ((μ^2)-(μ^2) Log[μ]+((μ^2)-(σ^2)) Log[(μ/(σ^2))]+2 (μ^2) Log[σ]+(σ^2) Log[Gamma[((μ^2)/(σ^2))]]+(-(μ^2)+(σ^2)) PolyGamma[0,((μ^2)/(σ^2))]) )//FullSimplify

(∂_σ) (-(μ^((μ^2)/(σ^2))) ((((μ/(σ^2))))^(-((μ^2)/(σ^2)))) (σ^(-2-((2 (μ^2))/(σ^2)))) ((μ^2)-(μ^2) Log[μ]+((μ^2)-(σ^2)) Log[(μ/(σ^2))]+2 (μ^2) Log[σ]+(σ^2) Log[Gamma[((μ^2)/(σ^2))]]+(-(μ^2)+(σ^2)) PolyGamma[0,((μ^2)/(σ^2))]) )//FullSimplify

Gamma in standard parameters mean μ dispersion κ

fn[x_,μ_,κ_]:=((((κ/μ)))^κ) (((E^(-x κ/μ)) (x^(-1+κ)) )/( Gamma[κ]))

TeXForm[((((κ/μ)))^κ) (((E^(-x κ/μ)) (x^(-1+κ)) )/( Gamma[κ]))]

Mean and variance

(∫_0^∞) ((((κ/μ)))^κ) (((E^(-x κ/μ)) (x^(-1+κ)) )/( Gamma[κ]))x \[DifferentialD]x//FullSimplify

(∫_0^∞) ((((κ/μ)))^κ) (((E^(-x κ/μ)) (x^(-1+κ)) )/( Gamma[κ]))(((x - μ))^2) \[DifferentialD]x //FullSimplify

Entropy

(∫_0^∞)-((((κ/μ)))^κ) (((E^(-x κ/μ)) (x^(-1+κ)) )/( Gamma[κ])) (κ Log[((κ/μ))] - Log[ Gamma[κ]]  - x (κ/μ) +(κ-1)Log[ x] )\[DifferentialD]x

-(-κ+Log[(κ/μ)]-Log[Gamma[κ]]+(-1+κ) PolyGamma[0,κ])//FullSimplify

κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ]//TeXForm

ESurfMuK=Plot3D[κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ], {μ,0.00,6},{κ,0.1,4},PlotRange->{-3,3},AxesLabel->{μ,κ},Mesh->None,LabelStyle->Directive[Medium,Bold]]

GamEntContMuK=ContourPlot[κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ], {μ,0.001,6},{κ,0.001,4},Contours->50,ColorFunction->(Hue[#]&), ClippingStyle->Automatic];

GradFlowMuK=GradientFieldPlot[κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ], {μ,0.001,6},{κ,0.001,4},GridLines->Automatic]

Show[GamEntContMuK,GradFlowMuK]

{(∂_μ)(κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ]),(∂_κ)(κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ])}//FullSimplify

{(1/μ),-(((-1+κ) (-1+κ PolyGamma[1,κ]))/κ)}/.{μ->x[t],κ->y[t]}

Solving grad flow curves in μ, κ

NDSolve[{x'[t]==(1/(x[t])),y'[t]==-(((-1+y[t]) (-1+PolyGamma[1,y[t]] y[t]))/(y[t])),x[0]==.1,y[0]==1},{x[t],y[t]},{t,0,50}]
k0110=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,50},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,6},{0,4}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==-(((-1+y[t]) (-1+PolyGamma[1,y[t]] y[t]))/(y[t])),x[0]==.1,y[0]==.1},{x[t],y[t]},{t,0,10}]
k0101=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,6},{0,4}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==-(((-1+y[t]) (-1+PolyGamma[1,y[t]] y[t]))/(y[t])),x[0]==.1,y[0]==2},{x[t],y[t]},{t,0,20}]
k0120=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,6},{0,4}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==-(((-1+y[t]) (-1+PolyGamma[1,y[t]] y[t]))/(y[t])),x[0]==.1,y[0]==3},{x[t],y[t]},{t,0,20}]
k0130=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,6},{0,4}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==-(((-1+y[t]) (-1+PolyGamma[1,y[t]] y[t]))/(y[t])),x[0]==1,y[0]==.1},{x[t],y[t]},{t,0,20}]
k1001=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,6},{0,4}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==-(((-1+y[t]) (-1+PolyGamma[1,y[t]] y[t]))/(y[t])),x[0]==2,y[0]==.1},{x[t],y[t]},{t,0,20}]
k2001=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,6},{0,4}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==-(((-1+y[t]) (-1+PolyGamma[1,y[t]] y[t]))/(y[t])),x[0]==3,y[0]==.1},{x[t],y[t]},{t,0,20}]
k3001=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,6},{0,4}}]

NDSolve[{x'[t]==(1/(x[t])),y'[t]==-(((-1+y[t]) (-1+PolyGamma[1,y[t]] y[t]))/(y[t])),x[0]==4,y[0]==.1},{x[t],y[t]},{t,0,20}]
k4001=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,20},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,6},{0,4}}]

GamEntFlowMuK=Show[k0110,k1001,k2001,k3001,k4001,k0101,k0130,k0101,k0120]

Show[GamEntContMuK,GradFlowMuK,GamEntFlowMuK]

GGradFlowMuK=GradientFieldPlot[κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ], {μ,0.001,1000},{κ,0.001,1000},GridLines->Automatic,AxesLabel->{μ,κ}];

GEntMuK=ContourPlot[κ-Log[(κ/μ)]+Log[Gamma[κ]]-(-1+κ) PolyGamma[0,κ], {μ,0.001,1000},{κ,0.001,1000},Contours->50,ColorFunction->(Hue[#]&), ClippingStyle->Automatic,GridLines->Automatic,AxesLabel->{{μ,κ}}];

Show[GEntMuK,GGradFlowMuK,AxesLabel->{μ,κ}];

McKay

Coords ((Subscript[α, 1]),(Subscript[σ, 12]),(Subscript[α, 2]))

m[x_,y_] :=((Subscript[α, 1])/(Subscript[σ, 12]))^((((Subscript[α, 1])+(Subscript[α, 2])))/2)  x^((Subscript[α, 1])-1) (y-x)^((Subscript[α, 2])-1) Exp[- (Sqrt[(Subscript[α, 1])/(Subscript[σ, 12])]) y]/(Gamma[(Subscript[α, 1])] Gamma[(Subscript[α, 2])])
m[x,y]

Integrate[x  m[x,y],{x,0,∞},{y,x,∞}]//PowerExpand//FullSimplify

Integrate[ y  m[x,y],{x,0,∞},{y,x,∞}]//PowerExpand//FullSimplify

Integrate[x y m[x,y],{x,0,∞},{y,x,∞}]//PowerExpand//FullSimplify

Entropy ((Subscript[α, 1]),(Subscript[σ, 12]),(Subscript[α, 2]))

Log[m[x,y]]m[x,y]//PowerExpand

Log[m[x,y]]m[x,y]//PowerExpand//FullSimplify

Integrate[Log[m[x,y]]m[x,y],{x,0,∞},{y,x,∞}]//PowerExpand//FullSimplify

ME[(Subscript[α, 1])_,(Subscript[σ, 12])_,(Subscript[α, 2])_]:=
(1/(Gamma[(Subscript[α, 1])] Gamma[(Subscript[α, 2])]))(-(Gamma[1+(Subscript[α, 1])] Gamma[(Subscript[α, 2])] (1+(((Log[(Subscript[α, 1])])/2)-((Log[(Subscript[σ, 12])])/2))-PolyGamma[0,(Subscript[α, 1])])+Gamma[(Subscript[α, 1])] (Gamma[1+(Subscript[α, 2])] (1+(((Log[(Subscript[α, 1])])/2)-((Log[(Subscript[σ, 12])])/2))-PolyGamma[0,(Subscript[α, 2])])+Gamma[(Subscript[α, 2])] (Log[Gamma[(Subscript[α, 1])]]-(-Log[Gamma[(Subscript[α, 2])]]+(1/2) ((Subscript[α, 1])+(Subscript[α, 2])) (Log[(Subscript[α, 1])]-Log[(Subscript[σ, 12])]))-2 (((Log[(Subscript[α, 1])])/2)-((Log[(Subscript[σ, 12])])/2))+PolyGamma[0,(Subscript[α, 1])]+PolyGamma[0,(Subscript[α, 2])]))) ((Subscript[α, 1]^(-((Subscript[α, 1])/2))) (Subscript[σ, 12]^(-(1/2) (-(Subscript[α, 1]))))) (Sqrt[(Subscript[σ, 12])]))//PowerExpand//FullSimplify

(1/(Gamma[(Subscript[α, 1])] Gamma[(Subscript[α, 2])]))(-(Gamma[1+(Subscript[α, 1])] Gamma[(Subscript[α, 2])] (1+(((Log[(Subscript[α, 1])])/2)-((Log[(Subscript[σ, 12])])/2))-PolyGamma[0,(Subscript[α, 1])])+Gamma[(Subscript[α, 1])] (Gamma[1+(Subscript[α, 2])] (1+(((Log[(Subscript[α, 1])])/2)-((Log[(Subscript[σ, 12])])/2))-PolyGamma[0,(Subscript[α, 2])])+Gamma[(Subscript[α, 2])] (Log[Gamma[(Subscript[α, 1])]]-(-Log[Gamma[(Subscript[α, 2])]]+(1/2) ((Subscript[α, 1])+(Subscript[α, 2])) (Log[(Subscript[α, 1])]-Log[(Subscript[σ, 12])]))-2 (((Log[(Subscript[α, 1])])/2)-((Log[(Subscript[σ, 12])])/2))+PolyGamma[0,(Subscript[α, 1])]+PolyGamma[0,(Subscript[α, 2])]))) ((Subscript[α, 1]^(-((Subscript[α, 1])/2))) (Subscript[σ, 12]^(-(1/2) (-(Subscript[α, 1]))))) (Sqrt[(Subscript[σ, 12])]))//PowerExpand//FullSimplify

(Subscript[α, 1]^(-((Subscript[α, 1])/2))) (((((2(ρ^2))/(1+(ρ^2)))))^((1+(Subscript[α, 1]))))(Log[((Subscript[α, 1])/((((((2(ρ^2))/(1+(ρ^2)))))^2)Gamma[(Subscript[α, 2])]Gamma[(Subscript[α, 1])]))]+PolyGamma[0,(Subscript[α, 1])] (-1+(Subscript[α, 1]))-(Subscript[α, 1])+PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))-(Subscript[α, 2])) //FullSimplify

(ρ^(-(ρ^2))) (((((2(ρ^2))/(1+(ρ^2)))))^((1+(ρ^2))))(Log[((ρ^2)/((((((2(ρ^2))/(1+(ρ^2)))))^2)Gamma[1-(ρ^2)]Gamma[(ρ^2)]))]+PolyGamma[0,(ρ^2)] (-1+(ρ^2))+
PolyGamma[0,1-(ρ^2)] (-(ρ^2))-1) //FullSimplify

TeXForm[%]

-(2^(1+1)) (((((ρ^4)/(((1+(ρ^2)))^2))))^((1/2) (1+1))) (1^(-(1/2))) (Log[((4 (ρ^4))/(((1+(ρ^2)))^2))]+Log[Gamma[1]]+Log[Gamma[(1/(ρ^2))-1]]-Log[1]-PolyGamma[0,1] (-1+1)+1-PolyGamma[0,(1/(ρ^2))-1] (-1+(1/(ρ^2))-1)+(1/(ρ^2))-1)//FullSimplify

TeXForm[((4 (ρ^4))/(((1+(ρ^2)))^2)) (Log[((4 (ρ^4)Gamma[(Subscript[α, 2])])/(((1+(ρ^2)))^2))]+(PolyGamma[0,(Subscript[α, 2])]-1) (-1+(Subscript[α, 2])))]

TeXForm[-((4 (ρ^4) (1+Log[((4 (ρ^4))/(((1+(ρ^2)))^2))]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2])))/(((1+(ρ^2)))^2))]

Plot3D[-((4 (ρ^4) (1+Log[((4 (ρ^4))/(((1+(ρ^2)))^2))]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2])))/(((1+(ρ^2)))^2)),{(Subscript[α, 2]),0.2,2},{ρ,0,0.99},PlotRange->All, AxesStyle->{9,9,9}]

D[-((4 (ρ^4) (1+Log[((4 (ρ^4))/(((1+(ρ^2)))^2))]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2])))/(((1+(ρ^2)))^2)),ρ]//FullSimplify

Plot3D[{0,((4 (ρ^4) (-1+PolyGamma[1,(Subscript[α, 2])] (-1+(Subscript[α, 2]))))/(((1+(ρ^2)))^2)),-(1/(((1+(ρ^2)))^3))16 (ρ^3) (2+Log[4]+Log[((ρ^4)/(((1+(ρ^2)))^2))]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2]))},{(Subscript[α, 2]),0,2},{ρ,0,0.99},PlotRange->{-7,5},PlotStyle->{RGBColor[0,1,0],RGBColor[1,0,0],RGBColor[0,0,1]}]

D[-((4 (ρ^4) (1+Log[((4 (ρ^4))/(((1+(ρ^2)))^2))]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2])))/(((1+(ρ^2)))^2)),(Subscript[α, 2])]//FullSimplify

Plot3D[{((4 (ρ^4) (-1+PolyGamma[1,(Subscript[α, 2])] (-1+(Subscript[α, 2]))))/(((1+(ρ^2)))^2))},{(Subscript[α, 2]),0,2},{ρ,0,0.99},PlotRange->{-7,5},PlotStyle->{RGBColor[1,0,0]}]

MCont=ContourPlot[-((4 (ρ^4) (1+Log[((4 (ρ^4))/(((1+(ρ^2)))^2))]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2])))/(((1+(ρ^2)))^2)),{(Subscript[α, 2]),0.01,2},{ρ,0.0,0.99},Contours->200,ColorFunction->(Hue[#]&),ClippingStyle->Automatic, GridLines->Automatic]

Needs["VectorFieldPlots`"]

MGrad=GradientFieldPlot[-((4 (ρ^4) (1+Log[((4 (ρ^4))/(((1+(ρ^2)))^2))]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2])))/(((1+(ρ^2)))^2)), {(Subscript[α, 2]), 0.01, 2}, {ρ, 0.01, .99},PlotRangeClipping->Automatic,Frame->True,GridLines->Automatic]

MContAndGrad=Show[MCont,MGrad]

Plot[-((4 (ρ^4) (1+Log[((4 (ρ^4))/(((1+(ρ^2)))^2))]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2])))/(((1+(ρ^2)))^2))/.(Subscript[α, 2])->0.4,{ρ,0.5,0.6},GridLines->Automatic]

V={{0.25,0.85},{0.3,0.66},{0.35,0.58},{0.4,0.53},{0.5,0.46},{0.6,0.425},{0.7,0.4},{0.8,0.38},{1,0.355},{1.1,0.345},{1.2,0.335},{1.3,0.33},{1.4,.325},{2,0.3}}

Clear[f]

f=Interpolation[V]

ME=Plot[{f[x],(1/((Sqrt[1+x])  )),(Sqrt[1-x])},{x,0.1,2},PlotRange->{{0,2.1},{0,0.97}},PlotStyle->{{AbsoluteThickness[7],RGBColor[0,0,0]},{AbsoluteThickness[2],RGBColor[0,0,0]},{AbsoluteThickness[2],RGBColor[0,0,0]}},Frame->True,GridLines->Automatic]

ML=ListPlot[V,PlotRange->{{0,2},{0,1}},PlotStyle->AbsolutePointSize[10],Frame->True,GridLines->Automatic]

Show[ME,ML]

NB To save 3D or colour graphics use Options, imagesize -> 600

Show[MContAndGrad,ME]

NDSolve[{x'[t]==((4 (-1+PolyGamma[1,x[t]] (-1+x[t])) ((y[t])^4))/(((1+((y[t])^2)))^2)),y'[t]==2+Log[4]+Log[Gamma[x[t]]]+Log[(((y[t])^4)/(((1+((y[t])^2)))^2))]-PolyGamma[0,x[t]] (-1+x[t])+x[t],x[0]==1.7,y[0]==.1},{x[t],y[t]},{t,0,5}]
m0405=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0.01,5},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,2},{0,1}}]

NDSolve[{x'[t]==((4 (-1+PolyGamma[1,x[t]] (-1+x[t])) ((y[t])^4))/(((1+((y[t])^2)))^2)),y'[t]==2+Log[4]+Log[Gamma[x[t]]]+Log[(((y[t])^4)/(((1+((y[t])^2)))^2))]-PolyGamma[0,x[t]] (-1+x[t])+x[t],x[0]==1,y[0]==.47},{x[t],y[t]},{t,0,5}]
m1005=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0.01,.405},Frame->True,GridLines->Automatic,PlotStyle->{Thickness[0.015],BoxRatios->{{1},{1}}},
PlotRange->{{0,2},{0,1}}]

NDSolve[{x'[t]==((4 (-1+PolyGamma[1,x[t]] (-1+x[t])) ((y[t])^4))/(((1+((y[t])^2)))^2)),y'[t]==2+Log[4]+Log[Gamma[x[t]]]+Log[(((y[t])^4)/(((1+((y[t])^2)))^2))]-PolyGamma[0,x[t]] (-1+x[t])+x[t],x[0]==1.5,y[0]==.3},{x[t],y[t]},{t,0,5}]
m1503=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0.01,.405},Frame->True,GridLines->Automatic,PlotStyle->{Thickness[0.015],BoxRatios->{{1},{1}}},
PlotRange->{{0,2},{0,1}}]

NDSolve[{x'[t]==((4 (-1+PolyGamma[1,x[t]] (-1+x[t])) ((y[t])^4))/(((1+((y[t])^2)))^2)),y'[t]==2+Log[4]+Log[Gamma[x[t]]]+Log[(((y[t])^4)/(((1+((y[t])^2)))^2))]-PolyGamma[0,x[t]] (-1+x[t])+x[t],x[0]==.5,y[0]==.5},{x[t],y[t]},{t,0,5}]
m0505=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0.01,.41},Frame->True,GridLines->Automatic,PlotStyle->{Thickness[0.015],BoxRatios->{{1},{1}}},
PlotRange->{{0,2},{0,1}}]

Coords ((Subscript[α, 1]),c,(Subscript[α, 2]))

m[x_,y_] :=c^((Subscript[α, 1])+(Subscript[α, 2]))  x^((Subscript[α, 1])-1) (y-x)^((Subscript[α, 2])-1) Exp[- c y]/(Gamma[(Subscript[α, 1])] Gamma[(Subscript[α, 2])])
m[x,y]

Integrate[x  m[x,y],{x,0,∞},{y,x,∞}]//PowerExpand//FullSimplify

Integrate[y  m[x,y],{x,0,∞},{y,x,∞}]//PowerExpand//FullSimplify

Integrate[x y  m[x,y],{x,0,∞},{y,x,∞}]//PowerExpand//FullSimplify

(((Subscript[α, 1]) (1+(Subscript[α, 1])+(Subscript[α, 2])))/(c^2))-( (1/(Gamma[(Subscript[α, 1])]))(c^(Subscript[α, 1]))(c^(-1-(Subscript[α, 1]))) Gamma[1+(Subscript[α, 1])])((((Subscript[α, 1])+(Subscript[α, 2]))/c) )//FullSimplify

Entropy ((Subscript[α, 1]),c,(Subscript[α, 2]))

GradientFieldPlot3D[(c^(-1-(Subscript[α, 1]))) (Sqrt[(Subscript[α, 1])]) (2 Log[c]-Log[Gamma[(Subscript[α, 1])]]-Log[Gamma[(Subscript[α, 2])]]+PolyGamma[0,(Subscript[α, 1])] (-1+(Subscript[α, 1]))-(Subscript[α, 1])+PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))-(Subscript[α, 2])),{(Subscript[α, 1]),.1,1}, {c,.1,1},{(Subscript[α, 2]),.1,1},VectorHeads->Automatic]

(1/(Gamma[(Subscript[α, 1])] Gamma[(Subscript[α, 2])]))(-(Gamma[1+(Subscript[α, 1])] Gamma[(Subscript[α, 2])] (1+(((Log[(Subscript[α, 1])])/2)-((Log[(Subscript[σ, 12])])/2))-PolyGamma[0,(Subscript[α, 1])])+Gamma[(Subscript[α, 1])] (Gamma[1+(Subscript[α, 2])] (1+(((Log[(Subscript[α, 1])])/2)-((Log[(Subscript[σ, 12])])/2))-PolyGamma[0,(Subscript[α, 2])])+Gamma[(Subscript[α, 2])] (Log[Gamma[(Subscript[α, 1])]]-(-Log[Gamma[(Subscript[α, 2])]]+(1/2) ((Subscript[α, 1])+(Subscript[α, 2])) (Log[(Subscript[α, 1])]-Log[(Subscript[σ, 12])]))-2 (((Log[(Subscript[α, 1])])/2)-((Log[(Subscript[σ, 12])])/2))+PolyGamma[0,(Subscript[α, 1])]+PolyGamma[0,(Subscript[α, 2])]))) ((Subscript[α, 1]^(-((Subscript[α, 1])/2))) (Subscript[σ, 12]^(-(1/2) (-(Subscript[α, 1]))))) (Sqrt[(Subscript[σ, 12])]))/.(Subscript[σ, 12])->((Subscript[α, 1])/(c^2))//FullSimplify

(c^(-1-(Subscript[α, 1]))) (Sqrt[(Subscript[α, 1])]) (2 Log[c]-Log[Gamma[(Subscript[α, 1])]]-Log[Gamma[(Subscript[α, 2])]]+PolyGamma[0,(Subscript[α, 1])] (-1+(Subscript[α, 1]))-(Subscript[α, 1])+PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))-(Subscript[α, 2]))/.(Subscript[α, 1])->1//FullSimplify

Plot3D[-((1-2 Log[c]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2]))/(c^2)), {c,0.1,5},{(Subscript[α, 2]),0.1,6}]

McFlow=GradientFieldPlot[-((1-2 Log[c]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2]))/(c^2)), {c,0.01,5},{(Subscript[α, 2]),0.01,5.8},PlotRangeClipping->Automatic,Frame->True,GridLines->Automatic]

McCont=ContourPlot[-((1-2 Log[c]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2]))/(c^2)), {c,0.01,5},{(Subscript[α, 2]),0.01,5.8},Contours->200,ColorFunction->(Hue[#]&),ClippingStyle->Automatic, GridLines->Automatic];

M1ContFlow=Show[McCont,McFlow]

Solving grad flow curves in (Subscript[M, 1]) where (Subscript[α, 1])=1

{(∂_c)(-((1-2 Log[c]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2]))/(c^2))),(∂_(Subscript[α, 2]))(-((1-2 Log[c]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2]))/(c^2)))}//FullSimplify

{((2 (2-2 Log[c]+Log[Gamma[(Subscript[α, 2])]]-PolyGamma[0,(Subscript[α, 2])] (-1+(Subscript[α, 2]))+(Subscript[α, 2])))/(c^3)),((-1+PolyGamma[1,(Subscript[α, 2])] (-1+(Subscript[α, 2])))/(c^2))}/.{c->x[t],(Subscript[α, 2])->y[t]}//FullSimplify

NDSolve[{x'[t]==((2 (2+Log[Gamma[y[t]]]-2 Log[x[t]]-PolyGamma[0,y[t]] (-1+y[t])+y[t]))/((x[t])^3)),y'[t]==((-1+PolyGamma[1,y[t]] (-1+y[t]))/((x[t])^2)),x[0]==.01,y[0]==1},{x[t],y[t]},{t,0,50}]
k00110=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,.69},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5},{0,5.8}}]

NDSolve[{x'[t]==((2 (2+Log[Gamma[y[t]]]-2 Log[x[t]]-PolyGamma[0,y[t]] (-1+y[t])+y[t]))/((x[t])^3)),y'[t]==((-1+PolyGamma[1,y[t]] (-1+y[t]))/((x[t])^2)),x[0]==.01,y[0]==2},{x[t],y[t]},{t,0,50}]
k00120=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,18.29},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5},{0,5.8}}]

NDSolve[{x'[t]==((2 (2+Log[Gamma[y[t]]]-2 Log[x[t]]-PolyGamma[0,y[t]] (-1+y[t])+y[t]))/((x[t])^3)),y'[t]==((-1+PolyGamma[1,y[t]] (-1+y[t]))/((x[t])^2)),x[0]==.01,y[0]==3},{x[t],y[t]},{t,0,97}]
k00130=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,110},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5},{0,5.8}}]

NDSolve[{x'[t]==((2 (2+Log[Gamma[y[t]]]-2 Log[x[t]]-PolyGamma[0,y[t]] (-1+y[t])+y[t]))/((x[t])^3)),y'[t]==((-1+PolyGamma[1,y[t]] (-1+y[t]))/((x[t])^2)),x[0]==.01,y[0]==4},{x[t],y[t]},{t,0,97}]
k00140=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,110},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5},{0,5.8}}]

NDSolve[{x'[t]==((2 (2+Log[Gamma[y[t]]]-2 Log[x[t]]-PolyGamma[0,y[t]] (-1+y[t])+y[t]))/((x[t])^3)),y'[t]==((-1+PolyGamma[1,y[t]] (-1+y[t]))/((x[t])^2)),x[0]==.01,y[0]==5},{x[t],y[t]},{t,0,97}]
k00150=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,110},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5},{0,5.8}}]

NDSolve[{x'[t]==((2 (2+Log[Gamma[y[t]]]-2 Log[x[t]]-PolyGamma[0,y[t]] (-1+y[t])+y[t]))/((x[t])^3)),y'[t]==((-1+PolyGamma[1,y[t]] (-1+y[t]))/((x[t])^2)),x[0]==.01,y[0]==1.5},{x[t],y[t]},{t,0,50}]
k00115=ParametricPlot[Evaluate[{x[t],y[t]}/.%],{t,0,5},Frame->True,GridLines->Automatic,PlotStyle->Thickness[0.015],PlotRange->{{0,5},{0,5.8}}]

M1Grad=Show[k00130,k00140,k00150,k00120,k00110,k00115]

Show[M1ContFlow,M1Grad]