The problem that i am facing is that of energy minimization. I have a function called totalenergy which is sum of two functions the elastic energy and the surface tension surfacetension.

totalenergy = k/2 (Integrate[r[angle], {angle, 0, Pi}] - c)^2 + 
Integrate[surfacetension[angle] r[angle]^2 Sin[angle], {angle, 0, Pi}]

k and c are constants in the above equation and angle is the angle in polar coordinates.

the user defines how surfacetension should should look like:

for example:

surfacetension[ang_] := PDF[NormalDistribution[0, 0.2], ang] + 0.5

Plot[surfacetension[x], {x, 0, Pi}, PlotRange -> {{0, Pi}, {0, Max[Pi]}}]

enter image description here

Now i wish to find the unknown profile r[angle] such that totalenergy is minimized for the known profile surfacetension[angle] for 0 <= angle <= Pi . This is precisely where I am stuck. Any help would be appreciated.


with marmots suggestion i tried the following code:

surfacetension[ang_ /; ang <= Pi] := PDF[NormalDistribution[0, 0.2], ang] + 0.5;
surfacetension[ang_ /; ang > Pi] := PDF[NormalDistribution[0, 0.2], 2 Pi - ang] + 0.5;
rsol[angle_, constant_] := r[angle] /. Solve[VariationalD[0.5 (r[angle] - 1)^2 +
surfacetension[angle] r[angle], r[angle], angle] == constant, r[angle]][[1]];
PolarPlot[rsol[\[Theta], 1], {\[Theta], 0, 2 Pi}]

enter image description here

however what i expect to get is the following

enter image description here

  • 1
    $\begingroup$ Are you looking for an exact solution? Otherwise you could discretise your unknown profile and you would then have a classical optimisation problem (which does not mean that it easy to solve; will probably depend on the number of unknowns), instead of a functional optimisation one. $\endgroup$ – anderstood Mar 28 '17 at 14:31
  • $\begingroup$ @anderstood actually right now i just need to crudely see what the function looks like. I am new to this area of Mathematica. Is it possible for you to post something as an answer $\endgroup$ – Ali Hashmi Mar 28 '17 at 14:34
  • 1
    $\begingroup$ Before being an "area of Mathematica", it is an area of mathematics. I suggest reading en.wikipedia.org/wiki/Euler%E2%80%93Lagrange_equation. I don't have time to investigate your problem right now but I might try later. $\endgroup$ – anderstood Mar 28 '17 at 14:41
  • 2
    $\begingroup$ @AliHashmi I agree that this is not a standard problem that can be solved immediately with EulerEquations. However, I believe that you can still compute the differential of the above energy functional. When I do this on paper, I get an equation telling me that r[alpha]=constant/(surfacetension[alpha Sin[alpha]) (no guarantee). Inserting this back and minimizing the result w.r.t. constant should give you the result (again no guarantee). $\endgroup$ – user46676 Mar 28 '17 at 16:57
  • 1
    $\begingroup$ Coding what I suggested might be more time than I have right now. In the unlikely event that I get some moments I'll give it a shot. $\endgroup$ – Daniel Lichtblau Mar 29 '17 at 15:58

I do not think that the following is to be considered an answer but may contain some steps towards an answer. The total energy is a functional of $r$, $E[r]$. You're looking for an extremum of this functional. The first term $E_1[r]$ is the square of a functional, $$ E_1[r]=\frac{k}{2}\left[\int\limits_0^\pi r(\alpha)\,\text{d}\alpha-c\right]^2=:\frac{k}{2}S[r]^2.$$ Using the chain rule, one obtains the differential $$ \frac{\delta E_1}{\delta r}(\alpha)=k\,S[r]\,\frac{\delta S}{\delta r}(\alpha)=k\,S[r]\,\cdot 1.$$ The second term can be treated with Mathematica's VariationalD. So the partial answer is the code

totalenergy := k/2 (NIntegrate[r[angle], {angle, 0, Pi}] - c)^2 + NIntegrate[surfacetension[angle] r[angle]^2 Sin[angle] , {angle, 0, Pi}];
dtotalenergy = k (Integrate[r[angle], {angle, 0, Pi}] - c) + VariationalD[surfacetension[angle] r[angle]^2 Sin[angle], r[angle], angle];
surfacetension[ang_] := PDF[NormalDistribution[0, 0.2], ang] + 0.5;
rsol[angle_, constant_] := r[angle] /.Solve[VariationalD[surfacetension[angle] r[angle]^2 Sin[angle], r[angle], 
   angle] == constant, r[angle]][[1]]; 

This seems to fix the functional form of $r$ to be
$$ r(\alpha) = \frac{0.5 \text{constant} \csc (\alpha)}{0.5\, +1.99471 e^{-12.5 \alpha^2}}.$$ The constant can be fixed by requiring that $\frac{\delta E}{\delta r}$ vanishes. However, I have problems with the integral. That's why I was proposing to redefine $r\to r\sqrt{\sin\alpha}$. As mentioned above, I'll be happy to erase this partial answer.

EDIT: replaced := before =. Also restarted the kernel, and the code went through.

  • $\begingroup$ thanks. can you quickly check the code for typos. because i can see := missing after rsol. can you run with the kernel restart and check if you get the answer that you have mentioned below $\endgroup$ – Ali Hashmi Mar 28 '17 at 21:55
  • $\begingroup$ what if instead of r[angle]^2 Sin[angle] we just have r[range]. I realized I might be mistaken with the inclusion of the square and the Sin term $\endgroup$ – Ali Hashmi Mar 29 '17 at 9:50
  • $\begingroup$ kindly see the edit. i made the surfacetension at angle = 0 and 2Pi the same $\endgroup$ – Ali Hashmi Mar 29 '17 at 11:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.