{NotebookFileName[],DateString[]}

style[
t_
]:=
Style[t,Bold,FontFamily->"Arial",FontSize->14]

$Assumptions={{x,y}∈Reals}

SeedRandom[61216281]

X (Y^2)

f1[
x_,
y_
]:=
x (y^2)

Maximize[{f1[x,y],x+y==200,x≥0,y≥0},{x,y}]

f1[
y_
]:=
f1[200-y,y]

f1[y]

D[f1[y],y]

D[f1[y],y]//Simplify

D[f1[x,y[x]],x]

D[f1[x,y[x]],x]/.y'[x]->-1

D[f1[x,y[x]],x]/.{y'[x]->-1,y[x]->y}

Solve[0==D[f1[x,y[x]],x]/.{y'[x]->-1,y[x]->y},{x,y}]

Block[{x=(200/3),y=(400/3)},
Graphics3D[{
Cuboid[{0,0,0},{x,y,y}],
White,
Text[style[X],{(x/2),0,0},{0,-1}],
Text[style[Y],{x,(y/2),y},{-1.5,0}],
Text[style[Y],{x,0,(y/2)},{1,0}],
Red,Thickness[0.05],
Line[{{x,(y/2),0},{x,(y/2),y},{0,(y/2),y},{0,(y/2),0},{x,(y/2),0}}]
}
]
]

Plot[f1[y],{y,0,200}
,PlotLegends->style/@{"X (Y^2), with X = 200-Y"}
,AxesLabel->style/@{"Y","V"}
]

(X^2)+2(Y^2)

f2[
x_,
y_
]:=
(x^2)+2(y^2)

Block[{xmin=-24,xmax=24,ymin=-24,ymax=24},
Show[
ContourPlot[{
f2[x,y]
},{x,xmin,xmax},{y,ymin,ymax}
,Contours->Automatic
,ContourLabels->True
,ColorFunction->(GrayLevel[1-#]&)
],
Plot[24-x,{x,xmin,xmax}
,PlotStyle->Directive[Thick,Red]
,PlotLegends->style/@{"X + Y = 24"}
,Frame->False
]
,PlotLabel->style["(X^2)+2(Y^2)"]
,FrameLabel->{{style["X"],None},{style["Y"],None}}
]
]

Block[{xmin=-24,xmax=24,ymin=-24,ymax=24,xs=16,ys=8},
Show[
ContourPlot[{
f2[x,y]==384
},{x,xmin,xmax},{y,ymin,ymax}
,Contours->Automatic
,ContourLabels->True
,PlotLegends->style/@{"(X^2) + 2(Y^2) = 384"}
,Epilog->{
PointSize[Large],
Point[{{16,8}}],
Text[style["("<>ToString[xs]<>", "<>ToString[ys]<>")"],{xs,ys},{1.2,0}]
}
],
Plot[24-x,{x,xmin,xmax}
,PlotStyle->Directive[Thick,Red]
,PlotLegends->style/@{"X + Y = 24"}
]
(*,PlotLabel->style["(X^2)+2(Y^2)"]*)
,FrameLabel->{{style["X"],None},{style["Y"],None}}
]
]

Minimize[{f2[x,y],x+y==24,x≥0,y≥0},{x,y}]

f2[24-y,y]//Expand

{f2[24,0],f2[0,24]}

D[f2[x,y[x]],x]

D[f2[x,y[x]],x]/.{y'[x]->-1,y[x]->y}

Solve[0==D[f2[x,y[x]],x]/.{y'[x]->-1,y[x]->y},{x,y}]

Plot[f2[24-y,y],{y,0,24}
,PlotRange->{0,All}
,PlotLegends->style/@{"X( ^2)+2(Y^2), with X = 24-Y"}
,AxesLabel->style/@{"Y",""}
]

Cylindrical can

$Assumptions=Union[{r>0,h>0},$Assumptions]

v3[
r_,
h_
]:=
π (r^2)h

a3[
r_,
h_
]:=
2π (r^2)+2π r h

Minimize[{a3[r,h],v3[r,h]==400,r>0,h>0},{r,h}]

Minimize[{a3[r,h],v3[r,h]==400,r>0,h>0},{r,h}]//N

D[v3[r,h[r]],r]

hprule=Solve[0==D[v3[r,h[r]],r],h'[r]][[1]]

D[a3[r,h[r]],r]

D[a3[r,h[r]],r]/.hprule/.h[r]->h

Solve[0==D[a3[r,h[r]],r]/.hprule/.h[r]->h,{r,h}]

rmin1=((400/(2π))^(1/3));
hmin1=2rmin1

a3[rmin1,hmin1]

a3[rmin1,hmin1]//N

Block[{r=3.99,h=7.98},
Graphics3D[{
Cylinder[{{0,0,0},{0,0,h}},r],
Thickness[0.02],,Red,
Line[{{0,0,h},{r,0,h}}],
Text[style["r"],{(r/2),0,h},{0,-1.5}],
Line[{{r,0,0},{r,0,h}}],
Text[style["h"],{r,0,(h/2)},{1.5,0}]

}
]
]

Solve[v3[r,h]==400,h]

v3[r,0]

a3[
r_
]:=Evaluate[
a3[r,h]/.Solve[v3[r,h]==400,h][[1]]
];
a3[r]

Plot[{
a3[r],
a3'[r]
},{r,0,10}
,PlotLegends->style/@{"A(r)","A'(r)"}
,AxesLabel->style/@{"r"}
,Epilog->{
PointSize[Large],
Point[{rmin1,a3[rmin1]}],
Text[style["r = "<>ToString[N@rmin1]<>" cm"],{rmin1,a3[rmin1]},{0,-4}],
Text[style["h = "<>ToString[N@2rmin1]<>" cm"],{rmin1,a3[rmin1]},{0,-2.5}],
Text[style["A = "<>ToString[N@a3[rmin1]]<>" (cm^2)"],{rmin1,a3[rmin1]},{0,-1}]
}
]

Disease outbreak

Number of patients as a function of time after disease outbreak.

p4[
t_
]:=
((t^2)/20)(E^(-t/3))

$Assumptions=Union[{t≥0},$Assumptions]

Maximize[{p4[t],0≤t≤20},t]

Fertilizer application

Yield on application of f kg of fertilizer. The problem here is just to maximize the reciprocal of a quadratic. This only makes sense if the quadratic is everywhere positive in the domain of the problem (which in this case is nonnegative numbers), since otherwise it will become unbounded in the vicinity of its roots. This being the case, it reduces to the problem of minimizing a quadratic. which is sort of too simple to be interesting.

y5[
f_
]:=
(2000/((f^2)-16f+100))

Maximize[{y5[f],f≥0},f]

$Assumptions=Union[{n>p>0},$Assumptions]

Sustainable harvest of deer

n is hundreds of deer. p is the breeding population in season i, n6(p) is the breeding population in the next year. The idea is to harvest n-p deer, leaving the population stable, and to pick p so as to maximize harvest.

n6[
p_
]:=
2Log[(p^2)+2p+1]

Maximize[{n6[p]-p,n>p>0},p]

D[n6[p]-p,p]

{p,n6[p]}/.Solve[0==D[n6[p]-p,p],p]

Plot[{n6[p],n6[p]-p,p},{p,0,5}
,PlotLegends->"Expressions"
]

Sustainable harvest of rabbits.

n7[
p_
]:=
16((Surd[p,4])^3)

Like previous example, but now p and n are number of individuals, I believe.

k7=Max[p/.Solve[n7[p]==p,p]]

Plot[{n7[p],n7[p]-p,p},{p,0,30000}
,PlotLegends->"Expressions"
]

D[n7[p]-p,p]

(12^4)

pmax1=p/.Solve[0==D[n7[p]-p,p],p][[1]]

{p,n7[p]-p,4*(12^3)}/.Solve[0==D[n7[p]-p,p],p]

Plot[n7[p]-p,{p,0,70000}
,AxesLabel->style/@{"N(t)","N'(t)"}
,PlotLegends->style/@{"N(t+1) - N(t)"}
,Epilog->{
PointSize[Large],
Point[{{k7,0}}],
Text[style["N = "<>ToString[k7]],{k7,0},{1.2,-1}]
}
]

Plot[Evaluate[{
n7[p]-p,
6000(n7'[p]-1)
}],{p,0,70000}
,PlotLegends->style/@{"H(N(t))","6000 H'(N(t))"}
,PlotRange->{-1000,7000}
,Epilog->{
PointSize[Large],
Point[{{pmax1,n7[pmax1]-pmax1}}],
Text[style["N(t) = "<>ToString[pmax1]],{pmax1,n7[pmax1]-pmax1},{0,1.5}],
Text[style["H(N(t)) = "<>ToString[n7[pmax1]-pmax1]],{pmax1,n7[pmax1]-pmax1},{0,3.5}]
}
]

Linear least-squares fitting

Assume a simple linear model, e.g. v(e)=(Subscript[k, cat])e for the velocity of an enzyme-catalyzed reaction as a function of added enzyme.

kcat1=100;
μ2=0.0;
σ2=100;
emax1=10;
n2=10

enzy2=Range[0,emax1,(emax1/n2)]

exact2=kcat1 enzy2

noise2=RandomReal[NormalDistribution[μ2,σ2],Length[exact2]]

data2=exact2+noise2

Round/@data2

ListPlot[{
Transpose[{enzy2,data2}]
}
,Joined->{False}
,PlotLegends->style/@{ "data"}
,AxesLabel->style/@{"[E]","v"}
,PlotStyle->ColorData[97][2]
,DataRange->{0,emax1}
,Filling->{2->{1}}
,FillingStyle->Directive[Thickness[0.01],Opacity[1]]
]

Block[{khats={50,110}},
ListPlot[{
Transpose[{enzy2,data2}],
Transpose[{enzy2,khats[[1]]enzy2}],
Transpose[{enzy2,khats[[2]]enzy2}]
}
,Joined->{False,True,True}
,PlotLegends->style/@{
"data",
ToString[khats[[1]]]<>"[E]",
ToString[khats[[2]]]<>"[E]"
}
,AxesLabel->style/@{"[E]","v"}
,PlotStyle->ColorData[97]/@{2,1,3}
,DataRange->{0,emax1}
(*,Filling->{1->{2},1->{3}}*)
,FillingStyle->Directive[Thickness[0.01],Opacity[1]]
]
]

ListPlot[{
Transpose[{enzy2,exact2}],
Transpose[{enzy2,data2}]
}
,Joined->{True, False}
,PlotLegends->{"exact linear", "data"}
,DataRange->{0,emax1}
,Filling->{2->{1}}
,FillingStyle->Directive[Thickness[0.01],Opacity[1]]
]

Φ((OverHat[k]))=(∑_(i=1)^n)((((OverHat[k]) (Subscript[e, i])-(Subscript[v, i])))^2)
Φ'((OverHat[k]))=(∑_(i=1)^n)2(Subscript[e, i])((OverHat[k])(Subscript[e, i])-(Subscript[v, i]))
=2(OverHat[k])(∑_(i=1)^n)(Subscript[e, i]^2)-2(∑_(i=1)^n)(Subscript[e, i])(Subscript[v, i])
(SuperStar[(OverHat[k])])=(((∑_(i=1)^n)(Subscript[e, i])(Subscript[v, i]))/((∑_(i=1)^n)(Subscript[e, i]^2)))

khat1=((Total[enzy2*data2])/(Total[(enzy2^2)]))

ListPlot[{
Transpose[{enzy2,khat1*enzy2}],
Transpose[{enzy2,data2}]
}
,Joined->{True, False}
,PlotLegends->{"least-squares linear fit", "data"}
,DataRange->{0,emax1}
,Filling->{2->{1}}
,FillingStyle->Directive[Thickness[0.01],Opacity[1]]
]

Nonlinear least-squares fitting

Generate a random exponential, y(t)=(e^(-λ t))+ϵ(t)

μ1=0.0;
σ1=0.05;
λ1=0.3;
n1=10

exact1=Table[
Exp[-λ1 t],
{t,0,n1}
]

noise1=RandomReal[NormalDistribution[μ1,σ1],Length[exact1]]

data1=exact1+noise1

ListPlot[{
data1
}
,Joined->{False}
,AxesLabel->style/@{"t (s)","[A]"}
,PlotLegends->style/@{ "[A]"}
,DataRange->{0,n1}
,Filling->{2->{1}}
,FillingStyle->Directive[Thickness[0.01],Opacity[1]]
]

ss1[
λ_
]:=Module[{curve,res},
curve=Table[Exp[-λ t],{t,0,n1}];
res=data1-curve;
res.res
]

D[(((Exp[-khat t[i]]-a[i]))^2),khat]

ss1[λ]

Plot[{
ss1[λ],
ss1'[λ]
},{λ,0,1}
,PlotLegends->{"SSR((OverHat[k]))","SSR'((OverHat[k]))"}
]

D[ss1[λ],λ]

D[ss1[λ],{λ,2}]//Expand

λhat1=λ/.FindRoot[D[ss1[λ],λ],{λ,0.1}]

Show[
ListPlot[{
data1
}
,Joined->{False}
,AxesLabel->style/@{"t (s)","[A]"}
,PlotLegends->style/@{ "[A]"}
,DataRange->{0,n1}
,Filling->{2->{1}}
,FillingStyle->Directive[Thickness[0.01],Opacity[1]]
],
Plot[Exp[-λhat1 t],{t,0,n1}
,PlotStyle->ColorData[97][2]
,PlotLegends->{"(E^(-(SuperStar[(OverHat[k])]) t))"}
]
]

Difficult case

f3[
x_
]:=
Cos[x]+Sin[x(Sqrt[2])]

Plot[{f3[x],f3'[x]},{x,-6π,6π}
,PlotLegends->"Expressions"
]

f3'[x]

Solve[0==f3'[x],x]

Entropy

s1[
p1_,
p2_,
p3_
]:=
-(p1 Log[p1]+p2 Log[p2]+p3 Log[p3]);
s1[p1,p2,p3]

s1[
p1_
]:=
s1[p1,1-p3s-p1,p3s];
s1[p1]

s1'[p1]

ss1[
p3_
]:=Module[{p12},
p12=((1-p3)/2);
-(2p12 Log[p12]+p3 Log[p3])
];ss1[p3]

Solve[0==ss1'[p3],p3]