Optimizing a Difference in Mathematica

Question

I have the following expressions:

$$ \frac{1}{2} \Bigg [1 + \sqrt{a^2 + (b e^{-\gamma(\tau)})^2 + (c e^{-\gamma(\tau)})^2} \Bigg ]; $$

f1=1/2 (1 + Sqrt[
    a^2 + (b Exp[-γ (\[tau)])^2 + (c Exp[-γ (\[tau)])^2])

and

$$ \frac{1}{2} \Bigg [1 + a^2 + e^{-\gamma(t)}(b^2 + c^2) \Bigg ]. $$

f2=1/2 (1 + a^2 + Exp[-γ (\[tau)] (b^2 + c^2))

I have already found these expressions such that for all values of $a, b$ and $c$ such that

$$ a, b, c \in [0, 1] \quad \text{such that} \quad |a|^2 + |b|^2 + |c|^2 = 1, $$

the first expression will always yield a higher value for all times, $\tau$.

Now, I am interested in figuring out certain values of $a, b$ and $c$ meeting the aforementioned condition for which the difference between the first and the second expressions is the largest. I have tried working with some sample numbers based on qualitative reasons.

Is there a specific method to go about doing this in Mathematica. Also, the difference will be different at different times due to the exponential term, which I'll get by running a numerical integration.

I'm not quite sure how to tackle this problem in Mathematica, if there's a specific method for it.

Edit:

Note that the values of the parameters (to follow) and the expression for gamma is given by the code:

G = 0.1;
\[Omega]c = 10;
\[Beta] = 1;

integralgamma[\[Omega]_, \[Tau]_] := 
4 G \[Omega] Exp[-\[Omega]/\[Omega]c] ((1 - 
   Cos[\[Omega] \[Tau]])/\[Omega]^(2)) Coth[\[Beta] \[Omega]/2];

I would run a numerical integration over omega first and then plot the graph with respect to tau, which is the time. The domain for the time is technically $(0, \infty)$ but for almost all cases, the curves of both expressions level off when tau (time) is in the range 10-20 seconds.

P.S: The centered LaTeX-ed expressions aren't written in Mathematica code. Read them as mathematical expressions.

Answer 1

Let us simplify your formulas manually first.

Let us denote $e^{-\gamma(t)}$ by e, and $a^2, b^2, c^2$ by a2, b2, c2.

Then

f1 = 1/2 (1 + Sqrt[a2 + b2 e^2 + c2 e^2]);

f2 = 1/2 (1 + a2 + e (b2 + c2));

You want to maximize f1 - f2 with the constraints $a^2+b^2+c^2=1$ and $a,b,c \in [0,1]$. From your notation, I also assume that $0 < e^{-\gamma(t)} < 1$.

result = Maximize[
  {f1 - f2, a2 + b2 + c2 == 1 && 0 < e < 1},
  {a2, b2, c2}
  ]

Refine[result, 0 < e < 1]

{(1 - 2 e + e^2)/(8 (1 + e)), 
 {a2 -> (1 + 3 e)/(4 (1 + e)), b2 -> 0, c2 -> (3 + e)/(4 (1 + e))}}

While I didn't give the constraints that $a,b,c \in [0,1]$, we can see that the result satisfies this.

---

If we use the method of Lagrange multipliers, and do the calculations with Mathematica, it turns out that the following is a solution for any 0 <= b2 <= 1/2:

{a2 -> (1 + 3 e)/(4 + 4 e), c2 -> (3 + e - 4 b2 (1 + e))/(4 (1 + e))}

Answer 2

If you introduce a spherical coordinates parametrization

$a=\sin(q)$, $b=\cos(q)\cos(\phi)$, $c=\cos(q)\sin(\phi)$,

you can write your problem in much simpler terms:

f1 = 1/2 (1 + Sqrt[Sin[q]^2 + k^2 Cos[q]^2])
f2 = 1/2 (1 + Sin[q]^2 + k Cos[q]^2)

where $k=\exp(\gamma(t))$

Now for a given value of $k=\exp(\gamma(t))$ you can do the 1d numerical maximization:

NMaximize[f1 - f2 /. {k -> 1}, {q}]
(*{0., {q -> -6.96782}}*)

Notice, $f1-f2$ do not depend on $b$ and $c$ in an independent fashion. They always enter in $b^2+c^2\equiv\cos^2(q)$ combination.

Then running

n = NMaximize[(f1 - f2) /. {k -> 1}, {q}][[2]]
o = FindInstance[((Cos[q] Cos[\[Phi]])^2 + (Cos[q] Sin[\[Phi]])^2 + 
      Sin[q]^2 == 1 && 0 < Cos[q] Cos[\[Phi]] < 1 && 
    0 < Cos[q] Sin[\[Phi]] < 1) /. n, \[Phi], Reals][[1]]

{ϕ -> -83.9}

Giving us

{a, b, c} = {Sin[q], Cos[q] Cos[\[Phi]], Cos[q] Sin[\[Phi]]} /. n /. o

{0.58074, 0.491245, 0.649168}

Answer 3

For a given a, b, and t and changing gamma, we can see that the diference between the two functions stops changing for large values of gamma.

  f1[a_, b_, c_, t_] := 
     1/2 (1 + Sqrt[
         a^2 + (b Exp[-\[Gamma] (t)])^2 + (c Exp[-\[Gamma] (t)])^2])
    f2[a_, b_, c_, t_] := 1/2 (1 + a^2 + Exp[-\[Gamma] (t)] (b^2 + c^2))
    v = {10., 15., 7.}
    {a, b, c} = v/Norm[v]
    a^2 + b^2 + c^2
    Plot[{f1[a, b, c, 2.], f2[a, b, c, 2.]}, {\[Gamma], 0, 10}]
    Clear[a, b, c]

So assuming a big value to gamma, we obtain:

Simplify[f1[a, b, c, 1.] - f2[a, b, c, 1.]] /. \[Gamma] -> 
   100000 // Chop

-(a^2/2)