Say we have some functional like the following: $H = (\partial_yf(y))^2 -w(y) f(y)^2 +f(y)^4/2$. This is the functional for the Gross Pitaevskii equation. Lets say $w(y)$, the trapping potential in this case, is some simple inverted Gaussian just so we have an explicit equation to work with: $w(y) = e^{-y^2}$.
My question is, what efficient methods are available in Mathematica for finding the function $f(y)$ that minimizes this functional? I will outline two methods that I have found that are problematic, but I would love to know some alternatives!
Variational Methods
We can use the variational methods. So I define my function and use NDSolve to solve the Euler-Lagrange equation. One issue I have here is that I don't know apriori what constraints I should use: since this gives me a 2nd order differential equation, I need 2 constraints. I set the derivative at y=0 to be 0, as I expect f(y) to be symmetric about y=0. Since I want the functional to be minimized, I use the fact that $f'(0)=0$ to also come up with the constraint that $f(0) = 1$. (To get this, I simple plugged the first constraint into the functional and solved for $f(0)$ that minimized the functional).
Needs["VariationalMethods`"]
w[y_]:=E^(-(y)^2)
sol = NDSolve[{EulerEquations[f'[y]^2 - w[y]*f[y]^2 + f[y]^4/2, f[y], y] , f'[0] == 0,f[0]==1}, f, {y, -5, 5}]
However, the result I am getting is certainly not a minimum of the functional. If one plots the above solution, you see that $f(y)$ blows up as |y| increases, which clearly can not be minimizing the functional because the functional has a term that goes at $f(y)^4$ which should keep $f(y)$ from blowing up. Am I doing something incorrectly here? I believe the issue is that it's trying to find an extrema rather than a Minima, which is what I am specifically looking for. However I don't know how to force it to minimize. If there were some way to simple use NMinimize to get a function, that would be ideal.
Discretize and Minimize
The other method I came up with may be the worst(most inefficient) way one could do this. The idea is the discretize y into pieces of size $\epsilon$, and treat each f[y] as an independent variable. One must also discretize the derivative as $f'[y] = (f[y]+f[y+\epsilon])/\epsilon$. Then you can directly minimize the functional as a function of N variables. The issue is that, if I want to know $f(y)$ with any precision, I need to make $\epsilon$ small, which means I have a very large number of variables. Thus, the minimization takes a very long time. Here is example code applying this to the above functional:
numbSite = 100;
var = Table[Subscript[\[Psi], i], {i, 0, numbSite - 1}];
\[Epsilon] = .1;
sites = Table[-\[Epsilon]*numbSite/2 + \[Epsilon]*i, {i, 1,numbSite}];
Ham2[var_?(VectorQ[#, NumericQ] &)] :=
Sum[((var[[i + 1]] - var[[i]])/\[Epsilon])^2, {i, 1, numbSite - 1}] -
Sum[(E^(-(-\[Epsilon]*numbSite/2 + \[Epsilon]*i)^2))*var[[i]]^2 , {i, 1, numbSite}] +
(1/2)*Sum[var[[i]]^4, {i, 1, numbSite}];
min = FindMinimum[{Ham2[var]},
Table[Subscript[\[Psi], i], {i, 0, numbSite - 1}],
MaxIterations -> 100];
p1 = ListPlot[
Table[{sites[[i]], var[[i]]^2 /. min[[2]]}, {i, 1, numbSite}]];
This gives closer to the result I am expecting. It's still very slow when I make the number of sites large, but at least it is something.
I would love to hear any suggestions you have on solving this problem! Functionals of this form are very common in physics (due to the usefulness of the Landau-Ginzburg theory) and it would be great to have more efficient ways of minimizing them.
