# 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}]][[1]][[
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. Commented Aug 30, 2020 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}]][[1]][[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);


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? Commented Aug 30, 2020 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. Commented Aug 30, 2020 at 4:32
• I see. That makes sense. Thanks a lot. Commented Aug 30, 2020 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. Commented Aug 30, 2020 at 0:09
• And how should one find the value of x where maxima occurs? Commented Aug 30, 2020 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. Commented Aug 30, 2020 at 0:42
• Any way of automating it? I want to do it for a bunch of different values. Commented Aug 30, 2020 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 )