How to find maxima of a numerical calculation

Question

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.

Bob Hanlon · Accepted Answer

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 *)

Hao Wang · Answer

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.

Hao Wang · Answer

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 )

How to find maxima of a numerical calculation

3 Answers

Add your own answers!

Ask a Question