2
$\begingroup$

I have a functional $S$,

$$S = \left( \int_{x_0}^{x_f} dx \frac{1}{z^d} \sqrt{-f(z,u) u'^2 - 2 u' z' +1} \right) + \frac{1}{z(x_0)^{d-1}}$$

where it is composed of two terms, i.e. an integral plus a constant term.

Varying $S$,

$$\delta S = \left( \int_{x_0}^{x_f} dx \rm{\;EOM} \right) + \frac{-(z' + f(z,u) u') \delta u - u' \delta z}{ z^d \sqrt{-f(z,u) u'^2 - 2 u' z' +1}}\Biggr|^{x_f}_{x_0} - \frac{d-1}{z(x_0)^d} \delta z(x_0)$$

The boundary term at $x_f$ vanishes since the points will be fixed. The only term left is at $x_0$ which is written as,

$$\frac{(z'(x_0) + f(z(x_0),u(x_0)) u'(x_0))}{ z(x_0)^d \sqrt{-f(z(x_0),u(x_0)) u'(x_0)^2 - 2 u'(x_0) z'(x_0) +1}} \delta u(x_0) + \left[\frac{u'(x_0)}{z(x_0)^d \sqrt{-f(z(x_0),u(x_0)) u'(x_0)^2 - 2 u'(x_0) z'(x_0) +1}} - \frac{d-1}{z(x_0)^d}\right]\delta z(x_0)$$

If we let $u(x_0)$ and $z(x_0)$ be movable boundary points then $\delta u(x_0) \neq 0$ and $\delta z(x_0) \neq 0$, so we have,

$$z'(x_0) + f(z(x_0),u(x_0)) u'(x_0) = 0$$

$$u'(x_0) - (d-1)\sqrt{-f(z(x_0),u(x_0)) u'(x_0)^2 - 2 u'(x_0) z'(x_0) +1} = 0$$

which actually could be combined to form,

$$u'(x_0) = \frac{-(d-1)^2 z'(x_0) \pm \sqrt{(d-1)^4 z'(x_0)^2 + 4(d-1)^2}}{2}$$

The boundary conditions of the problem are,

$$z(x_0) = z_s, z(x_f) = \epsilon = 10^{-5}$$

$$u'(x_0) = \frac{-(d-1)^2 z'(x_0) \pm \sqrt{(d-1)^4 z'(x_0)^2 + 4(d-1)^2}}{2}, u(x_f) = 1$$

Finding a solution to this problem will extremize $S$ for some chosen $z(x_0) = z_s$. The goal is to find the $z_s$ such that $S$ is minimized for all the set of solution that extremizes $S$.

The plot of $z(x)$ should look approximately as shown (ignore the values of the axes) where the blue dot is the location of $z(x_0) = z_s$,

Image

I have tried several ways to write a code, however, not only I had issues with the boundary value ODE but also I'm confused how and where to insert the minimization process involving $z(x_0)$. Any advice would be greatly appreciated.

ClearAll["Global`*"]
Needs["VariationalMethods`"]
d = 3;
u1 = 2;
m = (1/u[x] + 1/u1)^(d + 1);
f = 1 - m z[x]^(d + 1);
L = (Sqrt[-f u'[x]^2 - 2 u'[x] z'[x] + 1]/z[x]^d);
eulerlageq1 = EulerEquations[L, u[x], x];
eulerlageq2 = EulerEquations[L, z[x], x];
s = Solve[{eulerlageq1, eulerlageq2}, {u''[x], z''[x]}][[1]]//Simplify;
eq01 = {w'[x] == s[[1, 2]] /. u'[x] -> w[x], w[x] == u'[x]};
eq02 = {v'[x] == s[[2, 2]] /. z'[x] -> v[x], v[x] == z'[x]};

x0 = 10^-8;
xf = 10^-1;
eps = 10^-5; 
ztest[a_?NumericQ, b_?NumericQ, zs_?NumericQ] := First[z[xf] /. NDSolve[{eq01, eq02, z[x0] == zs, v[x0] == a, u[x0] == b, w[x0] == (-(d - 1)^2 a - Sqrt[(d - 1)^4 a^2 + 4 (d - 1)^2])/2}, {z, u}, {x, x0, xf}, Method -> "ExplicitEuler"]]
utest[a_?NumericQ, b_?NumericQ, zs_?NumericQ] := First[u[xf] /. NDSolve[{eq01, eq02, z[x0] == zs, v[x0] == a, u[x0] == b, w[x0] == (-(d - 1)^2 a - Sqrt[(d - 1)^4 a^2 + 4 (d - 1)^2])/2}, {z, u}, {x, x0, xf}, Method -> "ExplicitEuler"]]

sol = FindRoot[{ztest[a, b, zs] == eps, utest[a, b, zs] == 1}, {{a, -5}, {b, 10}}, WorkingPrecision -> 20]
S[zs_?NumericQ] := NIntegrate[Rationalize[L, 0] /. sol, {x, x0, xf}] + 1/zs^(d - 1)
FindMinimum[S[zs], {zs, 10}]
Improve this question
$\endgroup$
6
  • $\begingroup$ In your Mathematica code you write m = (1/u + 1/u1)^(d + 1), maybe you mean m = (1/u[x] + 1/u1)^(d + 1)? Since you say that you do not know how to deal this kind of problem, maybe it would be good to first consider a simpler problem that can be solved by hand, such as minimizing $S = u(1) + \int_0^1 ((u'(x))^2 + u(x)) dx$ with $u(0) = 0$. In fact, your question seems to be about math more than about Mathematica. Btw, if you can share the context where this problem shows up, that would be interesting. $\endgroup$
    – user293787
    Nov 24, 2022 at 17:10
  • $\begingroup$ @user293787 Thanks for the heads-up on the typo. Actually, I have considered simpler problems, i.e. has some conserved quantities so that the equation of motion can be written analytically. In that case, I could just take the derivative of $S$ with respect to $z_s$ and then find at which $z_s$ will $S$ be a minimum. In this case, however, there is no analytic form and I'd rather rely on FindMinimum to automate (simplify?) things. With regard to the context, I'm computing entanglement entropies where $S$ represents the classical area and the constant term is a quantum correction. $\endgroup$
    – mathemania
    Nov 24, 2022 at 18:15
  • $\begingroup$ @mathemania Could you map your region on (0,1) so that we can use wavelets as well? $\endgroup$
    – Alex Trounev
    Nov 25, 2022 at 4:58
  • $\begingroup$ Do you know any solution that satisfies the 3 boundary conditions for $z(x_f)$, $u(x_0)$, $u(x_f)$, no matter what $z(x_0)$? Btw, your eq01 still contains z' and your eq02 still contains u' (that is not necessarily a problem for NDSolve but I assume you intended to write the equation in standard 1st order form). Since $u,z>0$ I thought substituting $u = e^U$ and $z = e^Z$ could help, but maybe not. $\endgroup$
    – user293787
    Nov 25, 2022 at 8:42
  • $\begingroup$ @user293787 As of the moment I haven't solved any solution yet but if it can be solved I guess it should work for any $z(x_0)$ as long as it does allow the equations to diverge (refer to the equation above, you can see a $1/z^d$ factor). Regarding the two equations being solved using NDSolve, I'm having issues using the ShootingMethod as suggested in the MM post $(72725) (57262)$. $\endgroup$
    – mathemania
    Nov 25, 2022 at 15:03

1 Answer 1

Reset to default
4
$\begingroup$

First, we transform S to remove boundary value from optimization as follows $S = \left( \int_{x_0}^{x_f} dx \frac{1}{z^d} \sqrt{-f(z,u) u'^2 - 2 u' z' +1} \right) + \frac{1}{z(x_0)^{d-1}}$. We note, that $\frac{1}{z(x_0)^{d-1}}=\frac{1}{z(x_f)^{d-1}}+(d-1)\int_{x_0}^{x_f}\frac{z'dx}{z^d}$. Therefore, we should change Lagrange function as $L = \frac {1} {z^d}\sqrt {-f (z, u) u'^2 - 2 u' z' + 1} + \frac {(d - 1) z'} {z^d} $, then action is given by

$S = \int_ {x_ 0}^{x_f} Ldx + \frac {1} {\epsilon^{d - 1}} $

Last term is a constant, hence we can solve Euler equations to optimize action. We use the Euler wavelets collocation method

ClearAll["Global`*"]
Needs["VariationalMethods`"]
d = 3; zs = 10;
uc = 2; xf = 1/10; x0 = 1/xf 10^-8;
eps = 10^-5;
m = (1/u[x] + 1/uc)^(d + 1);
f = 1 - m ( z[x])^(d + 1);
L = (Sqrt[-(f u'[x]^2 + 2  u'[x] z'[x]) + 1]/( z[x])^
      d + (d - 1) z'[x]/z[x]^d);
eulerlageq1 = EulerEquations[L, u[x], x];
eulerlageq2 = EulerEquations[L, z[x], x];
s = Solve[{eulerlageq1, eulerlageq2}, {u''[x], z''[x]}][[1]] // 
   Simplify;
eq01 = u''[x] - s[[1, 2]] == 0;
eq02 = z''[x] - s[[2, 2]] == 0;
UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] := 
  Piecewise[{{2^(k/2) UE[m, 2^k t - 2 n + 1], (n - 1)/2^(k - 1) <= t <
       n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] := 
 Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 4; With[{k = k0, M = M0}, 
 nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; zcol = 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;
M = nn;
A = Array[a, {M}]; B = Array[b, M];
z2[x_] := A . Psi[x]; z1[x_] := A . int1[x] + a0; 
z0[x_] := A . int2[x] + a0 x + a1; u2[x_] := B . Psi[x]; 
u1[x_] := B . int1[x] + b0; u0[x_] := B . int2[x] + b0 x + b1;
var = Join[A, 
  B, {a0, a1, b0, b1}]; eqs = {u''[x] - s[[1, 2]], 
   z''[x] - s[[2, 2]]} /. {u''[x] -> u2[x]/xf^2, u'[x] -> u1[x]/xf, 
   u[x] -> u0[x], z''[x] -> z2[x]/xf^2, z'[x] -> z1[x]/xf, 
   z[x] -> z0[x]};
eq = Flatten[Table[eqs, {x, xcol}]];

bc = {z0[x0] == zs, z0[xs] == eps, u0[x0] == 15, 
    u0[xs] == 1} /. {xs -> 1};
action = 
 Table[L /. {u'[x] -> u1[x]/xf, u[x] -> u0[x], z'[x] -> z1[x]/xf, 
     z[x] -> z0[x]}, {x, xcol}] // Total; con = 
 Table[(-(f u'[x]^2 + 2 u'[x] z'[x]) + 1 >= 0) /. {u'[x] -> u1[x]/xf, 
    u[x] -> u0[x], z'[x] -> z1[x]/xf, z[x] -> z0[x]}, {x, xcol}];

sol2 = FindRoot[Join[Table[eq[[i]] == 0, {i, Length[eq]}], bc], 
  Table[{var[[i]], 1/10}, {i, Length[var]}], 
  MaxIterations -> Infinity, 
  Method -> {"Newton", "StepControl" -> "TrustRegion"}, 
  WorkingPrecision -> 30];

Visualization

{Plot[Evaluate[u0[x/xf] /. sol2], {x, x0, xf}, 
  AxesLabel -> {"x", "u"}, PlotRange -> All, PlotPoints -> 200], 
 Plot[Evaluate[z0[x/xf] /. sol2], {x, x0, xf}, 
  AxesLabel -> {"x", "z"}, PlotRange -> All, PlotPoints -> 200]}

Figure 1

Finally, we can evaluate action

action /. sol2

Out[]= 267.893419311596720291006268

To find minimum action with respect to zs we use Do loop as follows

Do[sol[j] = 
   FindRoot[
    Join[Table[eq[[i]] == 0, {i, Length[eq]}], 
     bc /. {zs -> 10 + (j - 5)/10}], 
    Table[{var[[i]], 1/10}, {i, Length[var]}], 
    MaxIterations -> Infinity, 
    Method -> {"Newton", "StepControl" -> "TrustRegion"}, 
    WorkingPrecision -> 30];, {j, 0, 10}]

Note that S is real in some region around zs=10, here we have table to look at 11 data with six real and 5 complex S (here we use unscaled S without constant part $1/\epsilon^2$)

S = 
 Table[{10 + (j - 5)/10, dx xf action /. sol[j]}, {j, 0, 10}]

Out[]= {{19/2, 
  1.57089633087946341494365634 + 
   0.000665399108477383332080644974 I}, {48/5, 
  1.57046931931274310853709700 + 
   0.000658066927631603946876114786 I}, {97/10, 
  1.57004300965468526183575503 + 
   0.000650226806636050815817267348 I}, {49/5, 
  1.56961754084009057627105968 + 
   0.000641969667063151355287558226 I}, {99/10, 
  1.67313437375896650260462628}, {10, 
  1.67433387069747950181878917}, {101/10, 
  1.67556219665128071840760756}, {51/5, 
  1.67682955598650942018621185}, {103/10, 
  1.67815118607052005825248523}, {52/5, 
  1.67955111804270498084313785}, {21/2, 
  1.52660724801448396934450389 + 0.0000609672885395141311715951931 I}}

Plot S shows that minimum is on the border between real and complex S at $9.8\le zs \le 9.9$

ListPlot[S]

Figure 2

Note, there is no big difference in 6 solutions as it shown in the picture below

{Plot[Evaluate[Table[u0[x/xf] /. sol[i], {i, 4, 9}]], {x, x0, xf}, 
  AxesLabel -> {"x", "u"}, PlotRange -> All], 
 Plot[Evaluate[Table[z0[x/xf] /. sol[i], {i, 4, 9}]], {x, x0, xf}, 
  AxesLabel -> {"x", "z"}, PlotRange -> All], 
 Plot[Evaluate[Table[z0[x/xf] /. sol[i], {i, 4, 9}]], {x, 0.035, .04},
   AxesLabel -> {"x", "z"}, PlotRange -> All, 
  PlotLegends -> Automatic]}

Figure 3

We can reproduce optimal solution with Haar wavelets as well. But in this case, there are no solutions with real action around zs=10. Nevertheless, plot for u looks like theoretical one, for example,

ClearAll["Global`*"]
Needs["VariationalMethods`"]
d = 3;
uc = 2; xf = 10^-1; x0 = 1/xf 10^-8; zs = 985/100;
eps = 10^-5;
m = (1/u[x] + 1/uc)^(d + 1);
f = 1 - m ( z[x])^(d + 1);
L = (Sqrt[-(f u'[x]^2 + 2  u'[x] z'[x]) + 1]/( z[x])^
      d + (d - 1) z'[x]/z[x]^d);
eulerlageq1 = EulerEquations[L, u[x], x];
eulerlageq2 = EulerEquations[L, z[x], x];
s = Solve[{eulerlageq1, eulerlageq2}, {u''[x], z''[x]}][[1]] // 
   Simplify;
eq01 = u''[x] - s[[1, 2]] == 0;
eq02 = z''[x] - s[[2, 2]] == 0;
J = 3; M = 2^J; dx = 1/(2*M); xl = Table[l*dx, {l, 0, 2*M}];
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2*M + 1}];
h1[x_] := Piecewise[{{1, 0 <= x <= 1}, {0, True}}];
p1[x_, n_] := (1/n!)*x^n; 
h[x_, k_, m_] := 
 Piecewise[{{1, 
    Inequality[k/m, LessEqual, x, Less, (1 + 2*k)/(2*m)]}, {-1, 
    Inequality[(1 + 2*k)/(2*m), LessEqual, x, Less, (1 + k)/m]}}, 0]
p[x_, k_, m_, n_] := 
 Piecewise[{{0, x < k/m}, {(-(k/m) + x)^n/n!, 
    Inequality[k/m, LessEqual, x, 
     Less, (1 + 2*k)/(2*m)]}, {((-(k/m) + x)^n - 
       2*(-((1 + 2*k)/(2*m)) + x)^n)/n!, (1 + 2*k)/(2*m) <= 
     x <= (1 + k)/
      m}, {((-(k/m) + x)^n + (-((1 + k)/m) + x)^n - 
       2*(-((1 + 2*k)/(2*m)) + x)^n)/n!, x > (1 + k)/m}}, 0]
var1 = Flatten[Table[a[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]];
var2 = Flatten[Table[b[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]];
z2[x_] := 
  Sum[a[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   a0*h1[x];
z1[x_] := 
  Sum[a[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   a0*p1[x, 1] + a1;
z0[x_] := 
  Sum[a[i, j]*p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   a0*p1[x, 2] + a1*x + a2;
u2[x_] := 
  Sum[b[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   b0*h1[x];
u1[x_] := 
  Sum[b[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   b0*p1[x, 1] + b1;
u0[x_] := 
  Sum[b[i, j]*p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   b0*p1[x, 2] + b1*x + b2;

var = Join[var1, var2, {a0, a1, a2, b0, b1, b2}];
eqs = {u''[x] == s[[1, 2]], 
    z''[x] == s[[2, 2]]} /. {u''[x] -> u2[x]/xf^2, u'[x] -> u1[x]/xf, 
    u[x] -> u0[x], z''[x] -> z2[x]/xf^2, z'[x] -> z1[x]/xf, 
    z[x] -> z0[x]};
eq = Flatten[Table[eqs, {x, xcol}]];

bc = {z0[x0] == zs, z0[xs] == eps, u0[x0] == 15, 
    u0[xs] == 1} /. {xs -> 1};
action = 
 Table[L /. {u'[x] -> u1[x]/xf, u[x] -> u0[x], z'[x] -> z1[x]/xf, 
     z[x] -> z0[x]}, {x, xcol}] // Total; con = 
 Table[(-(f u'[x]^2 + 2 u'[x] z'[x] zs)/xf^2 + 1 >= 0) /. {u'[x] -> 
     u1[x], u[x] -> u0[x], z'[x] -> z1[x], z[x] -> z0[x]}, {x, xcol}];

sol2 = FindRoot[Join[eq, bc], 
  Table[{var[[i]], 1/10}, {i, Length[var]}], 
  Method -> {"Newton", "StepControl" -> "TrustRegion"}, 
  MaxIterations -> Infinity, WorkingPrecision -> 30];

Visualization

{Plot[Evaluate[u0[x/xf] /. sol2], {x, x0, xf}, 
  AxesLabel -> {"x", "u"}, PlotRange -> All, PlotPoints -> 200], 
 Plot[Evaluate[z0[x/xf] /. sol2], {x, x0, xf}, 
  AxesLabel -> {"x", "z"}, PlotRange -> All, PlotPoints -> 200]}

Figure 4

Action has imaginary part

action /. sol2

Out[]= 340.6006157974879079324493482 + 0.109047774482325545497193142 I
Improve this answer
$\endgroup$
10
  • $\begingroup$ We know that the solution of the Euler-Lagrange equation extremizes the action $S$, my goal is to find which $z(x_0)$ will give a minimum out of all possible extremal solutions. You transformed the constant term $1/z(x_0)^{d-1}$ to $1/z(x_f)^{d-1}$ plus an integral is good since $1/z(x_f)^{d-1}$ now does not contribute to the variation, but how do you determine which $z(x_0)$ gives the minimum? In your code, you wrote zs=10. I think your code partly solves the problem where I don't have to worry about the term $1/z(x_0)^{d-1}$ since it's transformed into an integral. $\endgroup$
    – mathemania
    Nov 30, 2022 at 7:19
  • $\begingroup$ @mathemania I run last part of code, it takes a time to finish computation. The good news is that action at $z\le 9$, and $11\le z$ is complex, therefore the interval to find minimum is not so large. $\endgroup$
    – Alex Trounev
    Nov 30, 2022 at 7:36
  • $\begingroup$ The x0 you used is a bit different from mine. Also, in the line bc you used xs->1 instead of xf->10^-1. Why? I'm running the code but it's taking quite some time. $\endgroup$
    – mathemania
    Nov 30, 2022 at 8:42
  • $\begingroup$ @mathemania Since we use wavelets defined on $0\le x\le 1$ we map region $x0\le x\le xf$ to (0,1) using xf as a scale - see eqs definition. $\endgroup$
    – Alex Trounev
    Nov 30, 2022 at 9:00
  • $\begingroup$ I have found something a bit strange, when you evaluate for example, Evaluate[u0[x0/xf]/.sol2] you do not get u[x0]=15, what you get is u[x0] =14.980569299. I know this might have something to do with precision but I'm asking this since when I change my setup a little bit, using Haar wavelet produces a solution where it doesn't match the boundary condition by a considerable value. In this case, it just differs by an order of $10^{-2}$, but changing the setup can produce a difference of order $10^0$ which is bad. $\endgroup$
    – mathemania
    Dec 4, 2022 at 11:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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