# How to find maxima of a numerical calculation

Suppose I have the following code:

ClearAll["Global*"]
q := 1.6*10^-19; (* Electron charge in Coulomb *)
me := 9.1*10^-31; (* Free electron rest mass in kg *)
h := 6.63*10^-34;  (* Reduced Planck's constant in J.s *)
kb := 1.38*10^-23; (* Boltzmann constant in J/K *)
FD[d_, η_] := -PolyLog[
d + 1, -E^η];(* Defining the Fermi-Dirac integrals *)
Nc[d_, t_, α_, gs_, gv_, T_] :=
2*gs*gv*((kb*T)/α)^(d/t)* 1/(2*π^0.5)^d*1/t*Gamma[d/t]/
Gamma[d/2];  (* Effective band-edge DOS in d dimensions *)
n[d_, t_, α_, gs_, gv_, T_, ηF_] :=
Nc[d, t, α, gs, gv, T]*FD[(d - t)/t, ηF]*(100)^-d;
(* Effective SI carrier density in d dimensions in cm units*)
ηS[d_, t_, α_, gs_, gv_, T_, v_, nd_] :=
Quiet[Chop[
FindRoot[
1/2*(n[d, t, α, gs, gv, T, η] +
n[d, t, α, gs, gv, T, η - (q*v)/(kb*T)]) ==
nd, {η, 10000}]][][[
2]]]; (* Source Fermi Level at voltage v in d dimensions*)
J0[d_, t_, α_, gs_, gv_, T_] :=
gs*gv*q^2/h*((kb*T)/α)^((d - 1)/t)* 1/(2*π^0.5)^(d - 1)*(
kb*T)/q*Gamma[(d - 1 + t)/t]/Gamma[(d + 1)/2];
Jcore[d_, t_, α_, gs_, gv_, T_, v_, nd_] :=
FD[(d - 1)/t, ηS[d, t, α, gs, gv, T, v, nd]] -FD[(
d - 1)/t, ηS[d, t, α, gs, gv, T, v, nd] - (q*v)/(
kb*T)];
Jall[d_, t_, α_, gs_, gv_, T_, v_, nd_] :=
J0[d, t, α, gs, gv, T]*
Jcore[d, t, α, gs, gv, T, v, nd];
vinj[d_, t_, α_, gs_, gv_, T_, nd_, Eop_] :=
Jall[d, t, α, gs, gv, T, Eop/q, nd]/(q*nd);



On plotting vinj for some parameters mentioned below you should get the following plot:

v1 = LogLinearPlot[{vinj[2, 2, h^2/(4*π^2*2*0.2*me), 2, 1, 300,
n2d, 0.092*q]*10^-9}, {n2d, 10^10, 10^14}, Frame -> True,
FrameTicks -> Automatic,
FrameLabel -> {"Carrier Density (1/\!$$\*SuperscriptBox[\(cm$$, $$2\$$]\))", "\!$$\*SubscriptBox[\(v$$, $$inj$$]\) (\!$$\*SuperscriptBox[\ \(10$$, $$7$$]\)cm/s)"}, BaseStyle -> {FontSize -> 15},
PlotRange -> Full, AxesOrigin -> {Automatic, 0},
AspectRatio -> GoldenRatio, PlotStyle -> {Thick, Blue}] As you can see there is a maxima in the value of vinj as a function of n2d. How do I find this maxima numerically? Is there some function which can directly extract it? I tried FindMaxValue and NMaximize but both of them didn't work. Any suggestions are welcome.

• https://mathematica.stackexchange.com/a/26037/4999 might be helpful in getting FindMaxValue or FindMaximum to work. Aug 30 '20 at 2:59

ClearAll["Global*"]
q = 1.6*10^-19;(*Electron charge in Coulomb*)
me =
9.1*10^-31;(*Free electron rest mass in kg*)
h =
6.63*10^-34;(*Reduced Planck's constant in J.s*)
kb =
1.38*10^-23;(*Boltzmann constant in J/K*)

FD[d_, η_] := -PolyLog[
d + 1, -E^η];(*Defining the Fermi-Dirac integrals*)

Nc[d_, t_, α_, gs_, gv_, T_] :=
2*gs*gv*((kb*T)/α)^(d/t)*1/(2*π^(1/2))^d*1/t*
Gamma[d/t]/Gamma[d/2];(*Effective band-edge DOS in d dimensions*)
n[d_, t_, α_, gs_, gv_, T_, ηF_] :=
Nc[d, t, α, gs, gv, T]*FD[(d - t)/t, ηF]*(100)^-d;
(*Effective SI carrier density in d dimensions in cm units*)


Note that since ηS uses a numeric technique (FindRoot), its arguments should be restricted to numeric values using NumericQ

ηS[d_?NumericQ, t_?NumericQ, α_?NumericQ, gs_?NumericQ,
gv_?NumericQ, T_?NumericQ, v_?NumericQ, nd_?NumericQ] :=
Quiet[Chop[
FindRoot[1/
2*(n[d, t, α, gs, gv, T, η] +
n[d, t, α, gs, gv, T, η - (q*v)/(kb*T)]) == nd, {η,
10000}]][][]];(*Source Fermi Level at voltage v in d dimensions*)

J0[d_, t_, α_, gs_, gv_, T_] :=
gs*gv*q^2/h*((kb*T)/α)^((d - 1)/t)*1/(2*π^0.5)^(d - 1)*(kb*T)/q*
Gamma[(d - 1 + t)/t]/Gamma[(d + 1)/2];
Jcore[d_, t_, α_, gs_, gv_, T_, v_, nd_] :=
FD[(d - 1)/t, ηS[d, t, α, gs, gv, T, v, nd]] -
FD[(d - 1)/t, ηS[d, t, α, gs, gv, T, v, nd] - (q*v)/(kb*T)];
Jall[d_, t_, α_, gs_, gv_, T_, v_, nd_] :=
J0[d, t, α, gs, gv, T]*Jcore[d, t, α, gs, gv, T, v, nd];
vinj[d_, t_, α_, gs_, gv_, T_, nd_, Eop_] :=
Jall[d, t, α, gs, gv, T, Eop/q, nd]/(q*nd);


To find the maximum and argument use Maximize

Maximize[{vinj[2, 2, h^2/(4*π^2*2*(2/10)*me), 2, 1, 300,
n2d, (92/1000)*q]*10^-9, 10^10 < n2d < 10^14}, n2d]

(* {1.37927, {n2d -> 3.34063*10^12}} *)


For just the maximum use MaxValue

MaxValue[{vinj[2, 2, h^2/(4*π^2*2*(2/10)*me), 2, 1, 300,
n2d, (92/1000)*q]*10^-9, 10^10 < n2d < 10^14}, n2d]

(* 1.37927 *)


For just the argument of the maximum use ArgMax

ArgMax[{vinj[2, 2, h^2/(4*π^2*2*(2/10)*me), 2, 1, 300,
n2d, (92/1000)*q]*10^-9, 10^10 < n2d < 10^14}, n2d]

(* 3.34063*10^12 *)

• Thanks a lot. Though I don't understand what does NumericQ doing? Aug 30 '20 at 4:19
• NumericQ keeps the function from trying to evaluate with non-numeric arguments. See the "Applications" section of the documentation for NumericQ. Presumably, one of the reasons that you included Quiet in the function definition is that you got error messages when the function tried to evaluate with symbolic arguments. Aug 30 '20 at 4:32
• I see. That makes sense. Thanks a lot. Aug 30 '20 at 4:45

This is my solution:

vinjnew[x_] :=
Chop[vinj[2, 2, h^2/(4*\[Pi]^2*2*0.2*me), 2, 1, 300, x,
0.092*q]*10^-9]

vinjnew[#] & /@ Table[x, {x, 10^12, 10^13, 10^11}] // Max
vinjnew[#] & /@ Table[x, {x, 10^12, 10^13, 10^10}] // Max
vinjnew[#] & /@ Table[x, {x, 10^12, 10^13, 10^9}] // Max


The results of above codes are 1.37924, 1.37927, 1.37927 respectively. Hopefully you can try different precison like 10^8 to check the result.

I hope this can help.

• Thanks a lot. It should probably work. Aug 30 '20 at 0:09
• And how should one find the value of x where maxima occurs? Aug 30 '20 at 0:11
• I tried to find the x in a very primitive method, I think it is around 3.34* 10^12. While it is not very accurate. Aug 30 '20 at 0:42
• Any way of automating it? I want to do it for a bunch of different values. Aug 30 '20 at 0:45

I used this method to find the value of x

vinjnew2[x_] :=
Chop[vinj[2, 2, h^2/(4*\[Pi]^2*2*0.2*me), 2, 1, 300, x*10^12,
0.092*q]*10^-9];

Finder[x_] := Module[{para1, para2}, (
para1 = x;
If[Abs[vinjnew2[para1] - 1.37927] <= 0.000001, para2 = para1,
para2 = para1 + 0.01];
Return[para2])];

Nest[Finder, 3, 1000]


I think you may change the precison of x by changing the vinjnew function (I set x*10^12 inside the function) and the factor in the Finder (I set +0.01 per iterartion )