NDSolve with absolute value function

Question

I am trying to optimize a curve such that it minimalizes the following functional.

$$ E = \int_0^T dt \left| \left[ m x''(t) + C \sqrt{w^2 + 2 w x'(t) \sin{x(t)} + x'(t)^2} \left(w \sin{x(t)} + x'(t)\right)\right] x'(t) \right| $$

The variational derivative of this expression is a very long expression calculated by

eq = VariationalD[
  Abs[(m x''[t] + 
      c Sqrt[w^2 + 2 w x'[t] Sin[x[t]] + 
        x'[t]^2] (w  Sin[x[t]] + x'[t])) x'[t]], x[t], t]

The solution contains all kinds of derivatives of the Abs function, which I can imagine is tough to work with numerically. The NDSolve function encounters a 1/0 error.

pars = {m -> 80, c -> 0.25, w -> 2};
s = NDSolve[{(eq /. pars) == 0,
   x[0] == 0,
   x[1] == 2 \[Pi],
   x'[0] == 0,
   x''[0] == 0}, x, {t, 0, 1}]

I tried the following:

1. Substituting the absolute function with a differentiable absolute value function:

DifferentiableAbs[x_, h_] := Sqrt[x^2 + h^2]
eq = VariationalD[
  DifferentiableAbs[(m x''[t] + 
      c Sqrt[w^2 + 2 w x'[t] Sin[x[t]] + 
        x'[t]^2] (w  Sin[x[t]] + x'[t])) x'[t], h], x[t], t]

This takes a little while to calculate, but it does find a long expression in the end.

However, the NDSolve still gives a 1/0 error, which I had not expected. Any ideas how to avoid this?

2. Due to the underlying nature of the functional, I am quite certain that situations where the expression inside the absolute value in the integrand of the functional is negative can never be part of the solution where the functional is minimized. Therefore, I calculated the functional derivative without the absolute value, with the idea to try to find a way to restrict the integrand to be zero or positive when solving the equation.

eq = VariationalD[
  (m x''[t] + 
      c Sqrt[w^2 + 2 w x'[t] Sin[x[t]] + 
        x'[t]^2] (w  Sin[x[t]] + x'[t])) x'[t], x[t], t]

The functional derivative is quite elegant compared to the long expressions resulting from the derivative with the absolute value function.

NDSolve also finds a solution quite easily, but it contains situations where the integrand is negative.

pars = {m -> 80, c -> 0.25, w -> 2};
s = NDSolve[{(eq /. pars) == 0,
   x[0] == 0,
   x[1] == 2 \[Pi]}, x, {t, 0, 1}]
f[t_] := Evaluate[x[t] /. s /. pars]
P[t_] := ((m f''[t] + 
      c Sqrt[w^2 + 2 w f'[t] Sin[f[t]] + 
        f'[t]^2] (w Sin[f[t]] + f'[t])) f'[t] ) /. pars

I can't find a way to restrict the range of the solution in the solving process. Does anyone know if such an option exists? I did look into the WhenEvent functionality, but I couldn't get it to work.

Thanks in advance!

Update

Thanks for the help below! It is still not entirely solved though, unfortunately.

I updated the equation to give more realistic results.

$$ E = \int_0^T dt \left| \left[ m x''(t) + C \sqrt{w^2 + 2 w x'(t) \sin{\frac{2 \pi x(t)}{X}} + x'(t)^2} \left(w \sin{\frac{2 \pi x(t)}{X}} + x'(t)\right) + \frac{m g h'(x(t))}{\sqrt{1 + h'(x(t))^2}}\right] x'(t) \right| $$

with $x(0) = 0$ and $x(T) = X$.

Now I also choose more realistic parameters. I can find a solution to the non absolute action based equation for which the action is positive for all $t \in [0, T]$. I believe this means this is also a solution to the equation that is based on the absolute action.

h[x_, H_, X_] := H/2 (1 - Cos[(2 \[Pi] x)/X])
dhdx[x_, H_, X_] := H/X \[Pi] Sin[(2 \[Pi] x)/X]
L = (m x''[t] + 
     c Sqrt[w^2 + 2 w x'[t] Sin[(2 \[Pi] x[t])/X] + 
       x'[t]^2] (w Sin[(2 \[Pi] x[t])/X] + x'[t]) + 
     m g dhdx[x[t], H, X]/Sqrt[1 + dhdx[x[t], H, X]^2]) x'[t];
Labs = Abs[L];
pars = {m -> 80, g -> 9.8, c -> 1/4, w -> 4, X -> 30000, T -> 3600, 
   H -> 0, d -> 0};
eq = D[L, x[t]] - D[D[L, x'[t]], t] + D[D[L, x''[t]], t, t] /. pars //
    Simplify;

and

s = NDSolve[{eq == 0, x[0] == 0, x[T /. pars] == (X /. pars)}, 
  x[t], {t, 0, T /. pars}]
f[t_] := Evaluate[x[t] /. s[[1]]]
P[t_] := L  /. {x -> f} /. pars
ParametricPlot[{{f[t], f'[t]}, {f[t], X/T /. pars}}, {t, 0, 
  T /. pars}, AspectRatio -> 0.6, AxesLabel -> {"x", "x'"}, 
 PlotRange -> {{0, X /. pars}, {0, 12}}]
ParametricPlot[{f[t], P[t]}, {t, 0, T /. pars}, AspectRatio -> 0.6, 
 AxesLabel -> {"x", "P"}, PlotRange -> {{0, X /. pars}, {0, 250}}]

I compare the solution to the most trivial curve that satisfies the boundary conditions.

constant[t_] := X/T t;
NIntegrate[Labs /. {x -> constant} /. pars, {t, 0, T /. pars}]
NIntegrate[Labs /. {x -> f} /. pars, {t, 0, T /. pars}]

(*Out[441]= 609511.

Out[442]= 589477.*)

The solution NDSolve finds is better than the trivial curve.

Now, when I try to use the NMinimize method described below:

OEm[m_, x_] := 
 Sqrt[2 m + 
    1] Sum[(-1)^(m - k) x^k Binomial[m, k] Binomial[m + k, k], {k, 0, 
    m}]; UE[m_, t_] := OEm[m, t];
psi[k_, n_, m_, t_] := 
  Piecewise[{{2^((k - 1)/2) UE[m, 2^(k - 1) t - 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}]]]];
dt = 1/(nn); tl = Table[l*dt, {l, 0, nn}]; tcol = 
 Table[(tl[[l - 1]] + tl[[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;
A = Array[aa, nn];
x2[t_] := A . Psi[t];
x1[t_] := A . int1[t] + c1;
x0[t_] := A . int2[t] + c1 t + c2;


var = Join[{c1, c2}, A]; L0 = 
 Labs /. {x -> x0, x' -> x1, x'' -> x2}; action1 = 
 dt Sum[L0, {t, tcol}]; sol1 = 
 NMinimize[{action1 /. pars, {x0[0] == 0, 
    x0[T /. pars] == X /. pars}}, var, 
  Method -> "DifferentialEvolution"]
Plot[x0[t] /. sol1[[2]], {t, 0, T /. pars}]

This solution looks very linear. And if I compare its score against the NDSolve solution it leads me to believe that this is not a minimum. Or at least, I would have liked to find the NDSolve solution with this method as well.

NIntegrate[L0 /. sol1[[2]] /. pars, {t, 0, T /. pars}]

(*Out[474]= 609504.*)

Now, I believe the NDSolve solution is correct. It also makes sense from a physics viewpoint. However, I can tune $w$ and $h(x)$ in such a way that the action becomes negative, and then this is no longer a solution to the absoluted equation. Therefore, I would like to find a way to get this NMinimize method to work, or any other method.

Any ideas are very welcome!

Answer 1

We discussed a similar issue here. First, we should pay attention to the Lagrangian dependent on x''[t]. In this case equation of motion is given by

$\frac{\partial L}{\partial q} - \frac{d}{dt}\frac{\partial L}{\partial \dot q} + \frac{d^2}{dt^2}\frac{\partial L}{\partial \ddot q} = 0$

where $q=x(t)$. Using this equation we have

DifferentiableAbs[x_, h_] := Sqrt[x^2 + h^2];

L = DifferentiableAbs[(m x''[t] + 
     c Sqrt[w^2 + 2 w x'[t] Sin[x[t]] + x'[t]^2] (w Sin[x[t]] + 
        x'[t])) x'[t], h]; eq = 
 D[L, x[t]] - D[D[L, x'[t]], t] + D[D[L, x''[t]], t, t] /. {m -> 80, 
    c -> 1/4, w -> 2, h -> 0} // Simplify;

Solution

s = NDSolve[{eq == 0, x[0] == 0, x[1] == 2 \[Pi]}, x, {t, 0, 1}]

Visualization

Plot[x[t] /. s[[1]], {t, 0, 1}, Frame -> True, 
 FrameLabel -> {"t", "x"}]

Update 1. We can solve this problem using NMinimize and the Euler wavelets collocation method as follows

OEm[m_, x_] := 
 Sqrt[2 m + 
    1] Sum[(-1)^(m - k) x^k Binomial[m, k] Binomial[m + k, k], {k, 0, 
    m}]; UE[m_, t_] := OEm[m, t];
psi[k_, n_, m_, t_] := 
 Piecewise[{{2^((k - 1)/2) UE[m, 2^(k - 1) t - 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}]]]];
dt = 1/(nn); tl = Table[l*dt, {l, 0, nn}]; tcol = 
 Table[(tl[[l - 1]] + tl[[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;
A = Array[aa, nn];
x2[t_] := A . Psi[t];
x1[t_] := A . int1[t] + c1 ;
x0[t_] := A . int2[t] + c1 t + c2;


var = Join[{c1, c2}, A]; L0 = 
 L /. {x -> x0, x' -> x1, x'' -> x2}; action1 = 
 dt Sum[L0, {t, tcol}]; sol1 = 
 NMinimize[{action1 /. {m -> 80, c -> 1/4, w -> 2}, {x0[0] == 0, 
    x0[1] == 2 Pi}}, var, Method -> "DifferentialEvolution"]

Show[Plot[x[t] /. s[[1]], {t, 0, 1}, Frame -> True, 
  FrameLabel -> {"t", "x"}], 
 Plot[x0[t] /. sol1[[2]], {t, 0, 1}, PlotStyle -> {Red, Dashed}]]

Note that red dotted line computed with NMinimize considers with NDSolve solution around points t=0,1 only.

Update 2. Using numerical solutions s and sol we can compute L and $\int_0^1Ldt$ as folows

NIntegrate[L /. s[[1]] /. {m -> 80, c -> 1/4, w -> 2, h -> 0}, {t, 0, 1}]

(*Out[]= 1366.85*)
{Plot[L /. s[[1]] /. {m -> 80, c -> 1/4, w -> 2, h -> 0}, {t, 0, 1}, PlotLabel -> "NDSolve", AxesLabel -> {"t", "L"}], 
 Plot[L0 /. {m -> 80, c -> 1/4, w -> 2} /. sol[[2]], {t, 0, 1}, 
  PlotRange -> All, PlotLabel -> "NMinimize", 
  AxesLabel -> {"t", "L"}]}

As we can see from above solution computed with NDSolve is not minimum for action, but some maximum, while solution computed with NMInimize is minimum. Here the question arises, how to make the NDSolve calculate exactly the minimum solution?

Update 2 Second example with X=30000,T=3600 can be computed with the Euler wavelets collocation method as follows (we normalize x/X and t/T):

X = 30000; T = 3600;

h[x_, H_, X_] := H/2 (1 - Cos[(2 \[Pi] x)])
dhdx[x_, H_, X_] := H/X \[Pi] Sin[(2 \[Pi] x)]
L = (m x''[t]/X + 
     c Sqrt[w^2 + 2 w x'[t] Sin[(2 \[Pi] x[t])] + 
        x'[t]^2] (w Sin[(2 \[Pi] x[t])] + x'[t]) + 
     T^2/X^2 m g dhdx[x[t], H, X]/Sqrt[1 + dhdx[x[t], H, X]^2]) x'[
    t];
Labs = Abs[L];
pars = {m -> 80, g -> 9.8, c -> 1/4, w -> 4 T/X, H -> 0, d -> 0};

OEm[m_, x_] := 
 Sqrt[2 m + 
    1] Sum[(-1)^(m - k) x^k Binomial[m, k] Binomial[m + k, k], {k, 0, 
    m}]; UE[m_, t_] := OEm[m, t];
psi[k_, n_, m_, t_] := 
  Piecewise[{{2^((k - 1)/2) UE[m, 2^(k - 1) t - 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 = 5; M0 = 7;
With[{k = k0, M = M0}, 
  nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dt = 1/(nn); tl = Table[l*dt, {l, 0, nn}]; tcol = 
 Table[(tl[[l - 1]] + tl[[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;
A = Array[aa, nn];
x2[t_] := A . Psi[t];
x1[t_] := A . int1[t] + c1;
x0[t_] := A . int2[t] + c1 t + c2;


var = Join[{c1, c2}, A]; L0 = 
 Labs /. {x[t] -> x0[t], x'[t] -> x1[t], x''[t] -> x2[t]};
action1 = dt Sum[L0, {t, tcol}] /. pars; sol1 = 
 NMinimize[{action1, {x0[0] == 0, x0[1] == 1}}, var, 
  Method -> "DifferentialEvolution"]

Visualization together with linear solution 'X t/T`

Plot[{ X t/T, X x0[t/T] /. sol1[[2]]}, {t, 0, T}, 
 PlotStyle -> {{Automatic}, {Red, Dashed}}, AxesLabel -> Automatic]

Visualization together with NDSolve solution (blue solid line)