Finding the Minimum of a Four-Variable Function

Question

Greetings esteemed colleagues. I had previously posted regarding the minimisation of a function using Mathematica, and colleagues were kind enough to provide some useful suggestions that proved quite effective for my purposes. However, I have now encountered another issue that requires some undergraduate-level physics knowledge, though I shall attempt to keep matters straightforward or explain the relevant physics where needed.

Below is provided code that calculates the minimum reflectance of a two-layer system coated on a metal substrate. The angle of incidence (regarding the penetration of light through these systems) is also considered, but I have set this to zero for simplified calculations. The properties for the two-layer system are as follows:

First Layer:

• n1, for the refractive index.
• d1, for the thickness of layer 1.

Second layer:

• n2, for the refractive index.
• d2, for the thickness of layer 2.

(*This is a general code for S polarisation*)
ClearAll["Global`*"]
(*Primary initial Values*)
f = 390*10^9;(*Hz*)
\[Lambda] = 3*10^8/f;(*Metres*)
\[Sigma] = 2.5*10^6;(*Titanium*)
\[Epsilon] = 8.854*10^-12;
n = k = Sqrt[(\[Sigma])/(2*\[Epsilon]*2*\[Pi]*f)];
(*End of initial values*)

(*Setting up for calculation*)
(*Secondary initial value*)
\[Theta]1 = 0\[Degree];
\[Theta]2 = ArcSin[(n1*Sin[\[Theta]1])/n2];
\[Delta]1 = 2*\[Pi]*d1*n1*Cos[\[Theta]1]/\[Lambda];
\[Delta]2 = 2*\[Pi]*d2*n2*Cos[\[Theta]2]/\[Lambda];
\[Gamma]1 = n1*Cos[\[Theta]1];
\[Gamma]2 = n2*Cos[\[Theta]2];

(*Matrix calculation*)
matrixA = {{Cos[\[Delta]1],(I*Sin[\[Delta]1])/\[Gamma]1},{I*\[Gamma]1*Sin[\[Delta]1],Cos[\[Delta]1]}};
matrixB = {{Cos[\[Delta]2],(I*Sin[\[Delta]2])/\[Gamma]2},{I*\[Gamma]2*Sin[\[Delta]2],Cos[\[Delta]2]}};
matrixC = matrixA . matrixB;
Text["Here is Matrix BC:"];
matrixD = matrixC . {{1},{n-I*k}};
\[Beta] = matrixD[[1,1]];
\[Alpha] = matrixD[[2,1]];
Y = \[Alpha]/\[Beta];
Text["Here is r"];
r = (1-Y)/(1+Y);
R = Abs[r]^2;

(*Finding the minimum*)
Text["This is The minimum"]
NMinimize[{R, n1>0 && n2>0 && d1>0 && d2>0},{n1,d1,n2,d2}]

This is the result:

{0.941151, {n1 -> 2.69788, d1 -> 2.57474, n2 -> 0.991307, d2 -> 1.75953}}

Actually, the final line of code is of primary interest to me, as it calculates the minimum value of R with respect to n1, d1, n2, and d2. However, I have encountered an issue here. I have already attempted to manually manipulate the thicknesses of the layers (i.e. d1 and d2) and have also considered setting the thicknesses as λ/4 so that thickness becomes a function of the wavelength, which is an initial value. In this case, R is only a function of two variables (namely n1 and n2), and the result for this scenario is provided below:

{0.171573, {n1 -> 0.0101239, n2 -> 0.186525}}

In conclusion, I have found another "absolute minimum" lower than the previously recorded result; this indicates, as far as I can ascertain, that the initial code is unable to provide fully accurate results. I contemplated instructing Mathematica to consider n1 and d1 as one system, and n2 and d2 as another system, before determining the absolute minimum value. However, I am uncertain as to whether this approach would be successful.

Therefore, I would be most appreciative if colleagues could suggest any procedures or approaches (or a specific workflow) that may enable determining the absolute minimum value of R in a robust manner.

IMPORTANT UPDATE!

The numbers and minima are mathematically correct, but not physically! n1 and n2 must be greater than 1(I should have mentioned that in the condition for finding minimum but accidentally set that wrong, sorry). Also all of the variables have SI units, meaning that d1 and d2 are in metres which is physically wrong, it should be in the order of millimetres or micrometres. I've tried to set new conditions for NMinimize and here are the results:

NMinimize[{R, n1>1 && n2>1 && 0<d1<10^-4 && 0<d2<10^-4},{n1,n2,d1,d2}]

Result:

{0.128322, {n1 -> 1.24917, n2 -> 15.9991, d1 -> 0.0000583569, d2 -> 0.0000116479}}

Another condition:

NMinimize[{R, n1>1 && n2>1 && 0<d1<10^-3 && 0<d2<10^-3},{n1,n2,d1,d2}]

Result:

{0.901619, {n1 -> 3.5251, n2 -> 1.05889, d1 -> 0.000927416, d2 -> 0.000726193}}

Another condition:

NMinimize[{R, n1>1 && n2>1 && 0<d1<10^-4 && 0<d2<10^-3},{n1,n2,d1,d2}]

Result:

{0.789414, {n1 -> 1.02926, n2 -> 5.32453, d1 -> 2.91944*10^-6, d2 -> 0.000685964}}

To my knowledge, the results show that Mathematica cannot find the absolute minimum of the function R and it is bound to restricted conditions. I also assume that NMinimize is finding the local minima, not the absolute minimum. Can you help me with this?

Answer 1

NMinimize[{R, n1 > 0 && n2 > 0 && d1 > 0 && d2 > 0}, {n1, d1, n2, d2}, 
Method -> {"DifferentialEvolution", "ScalingFactor" -> 2.1},WorkingPrecision->15]

0.00728543143727867, {n1 -> 3.83485985429514, d1 -> 3.23665275424983, n2 -> 82.0912366363617, d2 -> 1.09616049856565}}

Addition. Omitting \[Theta]1 in your initial conditions, I obtain

NMinimize[{R,  n1 > 0 && n2 > 0 && d1 > 0 && d2 > 0 && \[Theta]1 >= 0 && 
\[Theta]1 <= 2*Pi}, {n1, d1, n2, d2, \[Theta]1},
 Method -> {"DifferentialEvolution", "ScalingFactor" -> 2.1},WorkingPrecision -> 15]

{0.00939550667791437, {n1 -> 25.0648945039521, d1 -> 4.85336449719848, n2 -> 206.930007948622, d2 -> 5.34646890841076, \[Theta]1 -> 6.19717877991621}}

Answer 2

Adding to @user64494's answer, another method that works quite well in this case is RandomSearch:

NMinimize[{R, n1 > 0 && n2 > 0 && d1 > 0 && d2 > 0}, {n1, d1, n2, d2},
  Method -> "RandomSearch", WorkingPrecision -> MachinePrecision]

{2.06244*10^-7, {n1 -> 1.38317, d1 -> 1.24808, n2 -> 28.9921, d2 -> 2.20152}}

And removing \[Theta1 =0 and adding it as a variable:

NMinimize[{R, 
  n1 > 0 && n2 > 0 && d1 > 0 && 
   d2 > 0 && \[Theta]1 >= 0 && \[Theta]1 <= 2*Pi}, {n1, d1, n2, 
  d2, \[Theta]1}, Method -> "RandomSearch", 
 WorkingPrecision -> MachinePrecision]

{8.46373*10^-6, {n1 -> 0.00425008, d1 -> 1.06997, n2 -> 1.19533, d2 -> 0.8329, \[Theta]1 -> 1.03731}}