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]}

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]

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]}

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]}

Action has imaginary part
action /. sol2
Out[]= 340.6006157974879079324493482 + 0.109047774482325545497193142 I
m = (1/u + 1/u1)^(d + 1)
, maybe you meanm = (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$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$eq01
still containsz'
and youreq02
still containsu'
(that is not necessarily a problem forNDSolve
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$NDSolve
, I'm having issues using theShootingMethod
as suggested in the MM post $(72725) (57262)$. $\endgroup$