(* These are automated verifications of the Bellman function properties. Some expressions are additionally transformed manually using the Mathematica commands. *)

Clear[f,g,p,q,r,t,s,gam,hess];
f[u_,v_,w_]=u^p gam[((v^q)/(u^p)),((w^r)/(u^p))];
MatrixForm[ReplaceAll[Expand[Simplify[(({{u^(-p/2+1), 0, 0}, {0, u^(-p/2+p/q), 0}, {0, 0, u^(-p/2+p/r)}})).ReplaceAll[Expand[(({{(∂_u)(∂_u)f[u,v,w], (∂_u)(∂_v)f[u,v,w], (∂_u)(∂_w)f[u,v,w]}, {(∂_u)(∂_v)f[u,v,w], (∂_v)(∂_v)f[u,v,w], (∂_v)(∂_w)f[u,v,w]-u}, {(∂_u)(∂_w)f[u,v,w], (∂_v)(∂_w)f[u,v,w]-u, (∂_w)(∂_w)f[u,v,w]}}))],
{v^q-> t u^p,w^r->s u^p,
v^(2q)->t^2 u^(2p),w^(2r)->s^2 u^(2p),
v^(-1+2q)->t^(2-1/q) u^(2p-p/q),w^(-1+2r)->s^(2-1/r) u^(2p-p/r),
v^(-1+q)->t^(1-1/q) u^(p-p/q),w^(-1+r)->s^(1-1/r) u^(p-p/r),
v^(-2+q)->t^(1-2/q) u^(p-2p/q),w^(-2+r)->s^(1-2/r) u^(p-2p/r),
v^(-2+2q)->t^(2-2/q) u^(2p-2p/q),w^(-2+2r)->s^(2-2/r) u^(2p-2p/r)}].(({{u^(-p/2+1), 0, 0}, {0, u^(-p/2+p/q), 0}, {0, 0, u^(-p/2+p/r)}}))]],
{((-1+q)/q)->1-1/q,((-2+q)/q)->1-2/q,((-1+r)/r)->1-1/r,((-2+r)/r)->1-2/r,-1+(1/q)+(1/r)->-1/p,p-(p/q)-(p/r)->1}]]

g[t_,s_]=(({{a+b t+c s}, {(a-(c/(p-1)))+b t+((c p)/(p-1)) s^(1-(1/p))}, {(a-((b+c)/(p-1)))+((b p)/(p-1)) t^(1-(1/p))+((c p)/(p-1)) s^(1-(1/p))}, {(a-((b+c)/(p-1)))+((b q)/2)t^((2/q)) s^((1/r)-(1/q))+(((c p)/(p-1))-((b p (q-r))/(2(p-1) r))) s^(1-(1/p))}, {(a-((b(q+r))/(2r(p-1))))+((b q^2)/(2 p(q-2))) t^((2/q))+ ((b q(q-r))/(2r(q-2))) t^((2/q)) s^(1-(2/q))+(c-((b(q-r))/(2r))) s}, {a+((b q)/(p(q-2))) t+((b q(q-r))/(2r(q-2))) t^((2/q)) s^(1-(2/q))+(c-((b(q-r))/(2r)))s}}));
h=(({{"1≤s≤t" }, {"s≤1≤t"}, {"s≤t≤1"}, {"t≤s≤1"}, {"t≤1≤s"}, {"1≤t≤s"}}));
Print[MatrixForm[Transpose[{g[t,s],h}]]];
regions={Graphics[{Thickness[.004],Line[{{0,0},{2,2}}]}],
Graphics[{Thickness[.004],Line[{{1,0},{1,2}}]}],
Graphics[{Thickness[.004],Line[{{0,1},{2,1}}]}],
Graphics[Text["t",{1.9,0.1}]],
Graphics[Text["s",{0.1,1.9}]],
Graphics[Text["1",{1.5,1.2}]],
Graphics[Text["2",{1.4,0.6}]],
Graphics[Text["3",{0.8,0.5}]],
Graphics[Text["4",{0.5,0.8}]],
Graphics[Text["5",{0.6,1.4}]],
Graphics[Text["6",{1.2,1.5}]]};
Show[regions,ImageSize->400,AspectRatio->Automatic,PlotRange->{{0,2},{0,2}},Axes->{True,True},Ticks->{None,None},TextStyle->{FontFamily->"Arial",FontSize->20}]

p=1/(1-1/q-1/r);
Print[Simplify[ReplaceAll[g[t,s][[1,1]],s->1]]==Simplify[ReplaceAll[g[t,s][[2,1]],s->1]]];
Print[Simplify[ReplaceAll[g[t,s][[2,1]],t->1]]==Simplify[ReplaceAll[g[t,s][[3,1]],t->1]]];
Print[Simplify[ReplaceAll[g[t,s][[3,1]],s->t]]==Simplify[ReplaceAll[g[t,s][[4,1]],s->t]]];
Print[Simplify[ReplaceAll[g[t,s][[4,1]],s->1]]==Simplify[ReplaceAll[g[t,s][[5,1]],s->1]]];
Print[Together[Simplify[ReplaceAll[g[t,s][[5,1]],t->1]]]==Together[Simplify[ReplaceAll[g[t,s][[6,1]],t->1]]]];
Print[Simplify[ReplaceAll[g[t,s][[6,1]],s->t]]==Simplify[ReplaceAll[g[t,s][[1,1]],s->t]]];

Print[Simplify[ReplaceAll[(∂_t)g[t,s][[1,1]],s->1]]==Simplify[ReplaceAll[(∂_t)g[t,s][[2,1]],s->1]]];
Print[Simplify[ReplaceAll[(∂_t)g[t,s][[2,1]],t->1]]==Simplify[ReplaceAll[(∂_t)g[t,s][[3,1]],t->1]]];
Print[Simplify[ReplaceAll[(∂_t)g[t,s][[3,1]],s->t]]==Simplify[ReplaceAll[(∂_t)g[t,s][[4,1]],s->t]]];
Print[Simplify[ReplaceAll[(∂_t)g[t,s][[4,1]],s->1]]==Simplify[ReplaceAll[(∂_t)g[t,s][[5,1]],s->1]]];
Print[Simplify[ReplaceAll[(∂_t)g[t,s][[5,1]],t->1]]==Simplify[ReplaceAll[(∂_t)g[t,s][[6,1]],t->1]]];
Print[Simplify[ReplaceAll[(∂_t)g[t,s][[6,1]],s->t]]==Simplify[ReplaceAll[(∂_t)g[t,s][[1,1]],s->t]]];

Print[Simplify[ReplaceAll[(∂_s)g[t,s][[1,1]],s->1]]==Simplify[ReplaceAll[(∂_s)g[t,s][[2,1]],s->1]]];
Print[Simplify[ReplaceAll[(∂_s)g[t,s][[2,1]],t->1]]==Simplify[ReplaceAll[(∂_s)g[t,s][[3,1]],t->1]]];
Print[Simplify[ReplaceAll[(∂_s)g[t,s][[3,1]],s->t]]==Simplify[ReplaceAll[(∂_s)g[t,s][[4,1]],s->t]]];
Print[Simplify[ReplaceAll[(∂_s)g[t,s][[4,1]],s->1]]==Simplify[ReplaceAll[(∂_s)g[t,s][[5,1]],s->1]]];
Print[Simplify[ReplaceAll[(∂_s)g[t,s][[5,1]],t->1]]==Simplify[ReplaceAll[(∂_s)g[t,s][[6,1]],t->1]]];
Print[Simplify[ReplaceAll[(∂_s)g[t,s][[6,1]],s->t]]==Simplify[ReplaceAll[(∂_s)g[t,s][[1,1]],s->t]]];

Clear[t,s,p,q,r];
For[j=6,j≤6,j++,
gam[t_,s_]=g[t,s][[j,1]];
Print["Region ",j,": ",h[[j,1]]];
Print["Minor 1x1:"];
hess[t_,s_]=(-p gam[t,s]+(p^2) gam[t,s]+p s Derivative[0,1][gam][t,s]-(p^2) s Derivative[0,1][gam][t,s]+(p^2) (s^2) Derivative[0,2][gam][t,s]+p t Derivative[1,0][gam][t,s]-(p^2) t Derivative[1,0][gam][t,s]+2 (p^2) s t Derivative[1,1][gam][t,s]+(p^2) (t^2) Derivative[2,0][gam][t,s]);
minexp=Expand[FullSimplify[ReplaceAll[Expand[Simplify[hess[t,s]]],{1/p->1-1/q-1/r}]]];
Print[minexp];
Print["Minor 2x2:"];
hess[t_,s_]=(({{-p gam[t,s]+(p^2) gam[t,s]+p s Derivative[0,1][gam][t,s]-(p^2) s Derivative[0,1][gam][t,s]+(p^2) (s^2) Derivative[0,2][gam][t,s]+p t Derivative[1,0][gam][t,s]-(p^2) t Derivative[1,0][gam][t,s]+2 (p^2) s t Derivative[1,1][gam][t,s]+(p^2) (t^2) Derivative[2,0][gam][t,s], -p q s (t^(1-(1/q))) Derivative[1,1][gam][t,s]-p q (t^(2-(1/q))) Derivative[2,0][gam][t,s]}, {-p q s (t^(1-(1/q))) Derivative[1,1][gam][t,s]-p q (t^(2-(1/q))) Derivative[2,0][gam][t,s], -q (t^(1-(2/q))) Derivative[1,0][gam][t,s]+(q^2) (t^(1-(2/q))) Derivative[1,0][gam][t,s]+(q^2) (t^(2-(2/q))) Derivative[2,0][gam][t,s]}}));
minexp=Expand[FullSimplify[ReplaceAll[Expand[Simplify[Det[hess[t,s]]]],{1/p->1-1/q-1/r}]]];
Print[minexp];
Print["Determinant with -:"];
hess[t_,s_]=(({{-p gam[t,s]+(p^2) gam[t,s]+p s Derivative[0,1][gam][t,s]-(p^2) s Derivative[0,1][gam][t,s]+(p^2) (s^2) Derivative[0,2][gam][t,s]+p t Derivative[1,0][gam][t,s]-(p^2) t Derivative[1,0][gam][t,s]+2 (p^2) s t Derivative[1,1][gam][t,s]+(p^2) (t^2) Derivative[2,0][gam][t,s], -p q s (t^(1-(1/q))) Derivative[1,1][gam][t,s]-p q (t^(2-(1/q))) Derivative[2,0][gam][t,s], -p r (s^(2-(1/r))) Derivative[0,2][gam][t,s]-p r (s^(1-(1/r))) t Derivative[1,1][gam][t,s]}, {-p q s (t^(1-(1/q))) Derivative[1,1][gam][t,s]-p q (t^(2-(1/q))) Derivative[2,0][gam][t,s], -q (t^(1-(2/q))) Derivative[1,0][gam][t,s]+(q^2) (t^(1-(2/q))) Derivative[1,0][gam][t,s]+(q^2) (t^(2-(2/q))) Derivative[2,0][gam][t,s], -1+q r (s^(1-(1/r))) (t^(1-(1/q))) Derivative[1,1][gam][t,s]}, {-p r (s^(2-(1/r))) Derivative[0,2][gam][t,s]-p r (s^(1-(1/r))) t Derivative[1,1][gam][t,s], -1+q r (s^(1-(1/r))) (t^(1-(1/q))) Derivative[1,1][gam][t,s], -r (s^(1-(2/r))) Derivative[0,1][gam][t,s]+(r^2) (s^(1-(2/r))) Derivative[0,1][gam][t,s]+(r^2) (s^(2-(2/r))) Derivative[0,2][gam][t,s]}}));
minexp=Expand[FullSimplify[ReplaceAll[Expand[Simplify[Det[hess[t,s]]]],{1/p->1-1/q-1/r}]]];
Print[minexp];
Print["Determinant with +:"];
hess[t_,s_]=(({{-p gam[t,s]+(p^2) gam[t,s]+p s Derivative[0,1][gam][t,s]-(p^2) s Derivative[0,1][gam][t,s]+(p^2) (s^2) Derivative[0,2][gam][t,s]+p t Derivative[1,0][gam][t,s]-(p^2) t Derivative[1,0][gam][t,s]+2 (p^2) s t Derivative[1,1][gam][t,s]+(p^2) (t^2) Derivative[2,0][gam][t,s], -p q s (t^(1-(1/q))) Derivative[1,1][gam][t,s]-p q (t^(2-(1/q))) Derivative[2,0][gam][t,s], -p r (s^(2-(1/r))) Derivative[0,2][gam][t,s]-p r (s^(1-(1/r))) t Derivative[1,1][gam][t,s]}, {-p q s (t^(1-(1/q))) Derivative[1,1][gam][t,s]-p q (t^(2-(1/q))) Derivative[2,0][gam][t,s], -q (t^(1-(2/q))) Derivative[1,0][gam][t,s]+(q^2) (t^(1-(2/q))) Derivative[1,0][gam][t,s]+(q^2) (t^(2-(2/q))) Derivative[2,0][gam][t,s], 1+q r (s^(1-(1/r))) (t^(1-(1/q))) Derivative[1,1][gam][t,s]}, {-p r (s^(2-(1/r))) Derivative[0,2][gam][t,s]-p r (s^(1-(1/r))) t Derivative[1,1][gam][t,s], 1+q r (s^(1-(1/r))) (t^(1-(1/q))) Derivative[1,1][gam][t,s], -r (s^(1-(2/r))) Derivative[0,1][gam][t,s]+(r^2) (s^(1-(2/r))) Derivative[0,1][gam][t,s]+(r^2) (s^(2-(2/r))) Derivative[0,2][gam][t,s]}}));
minexp=Expand[FullSimplify[ReplaceAll[Expand[Simplify[Det[hess[t,s]]]],{1/p->1-1/q-1/r}]]];
Print[minexp];
];

(* We need to make each expression ≥0. This is achieved by
taking b arbitrary
taking c large enough (depending on b,p,q,r)
taking a large enough (depending on b,c,p,q,r). *)
(* The dominant terms (regarding the exponents of s,t) are shaded in blue. In each case and in each expression there is always only one dominant term. *)
(* We need to take care that the factor of a in the dominant term can be chosen strictly larger than the one in the term we want to dominate. This is not the case for the terms shaded in green, but these terms can be made ≥0, so they need not be dominated. *)