# NDSolve with absolute value function

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?

1. 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.

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!

• I believe that the Calculus of Variations can be applied only to functionals that are continuously differentiable with respect to the function and its derivatives. So, only your second approach is valid. Also, the term m x''[t] x'[t] drops out of the calculation. As a result, m and c drop out too. Is this your intent? Aug 13 at 0:57
• Why not to use NMinimize? Aug 13 at 1:41
• @bbgodfrey Yes, I know m drops out (c doesn't, though). The integrand is actually the physical quantity power (force times velocity). For that reason I left it in. Aug 13 at 11:29
• @AlexTrounev Is NMinimize able to find a curve? It looks like it only finds single values. Aug 13 at 11:31
• @JanMulder Yes, we can find a curve using NMinimize. See Update 1 to my answer. Aug 13 at 12:50

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)

• The first solution seems to work on first glance. However, the final solution is the same as the solution without the absolute sign. I think the NDSolve rewrites all the Sqrt[x^2] as x, therefore not really applying the absolute functions. Aug 13 at 15:15
• The problem with variational principle is that action is not really has a minimum on the solution of eq==0 but it is rather stationary point. If we compute NIntegrate[L /. s[[1]] /. {m -> 80, c -> 1/4, w -> 2, h->0}, {t, 0, 1}], then we have 1366.85, while with NMinimize` we have about 0.0000196266. The last value probably is minimum, but 1366.85 looks like a maximum. Aug 13 at 15:42
• Thanks for your efforts to solve this, you already helped me a lot. Much appreciated! I believe the NDSolve solution is the solution to the non absolute equation. It looks like only the action is absoluted. The NMinimize solution also doesn't look like it's correct. I provided an update of the problem in the initial description. Aug 20 at 9:05