# NDSolve problem with conservation laws

I am not an experienced Mathematica user, but decided to use it this time because of the in-buit special functions that it has.

My problem is a solution of the following ODE

y''[x]=-f(y[x]),  y'[0]=y'[L]=c,


with $$L>0$$, $$c>0$$, and periodic rhs $$f(y+2\pi)=f(y)=-\Im\{\text{Li}_{2}(-e^{iy})\}$$. I tried to solve the equation with the in-built NDSolve function and encountered the following problem.

Using the first integral of the above equation we find that

0.5*(y'[x])^2+F(y[x])=C,


where $$C$$ is a constant, and $$F(y)=F(y+2\pi)$$ such that $$F'(y)=f(y)$$. Using the boundary conditions we infer that

F(y[0])=F(y[L]),


implying that y[L]=y[0]+2*Pi*n. While solving the equation numerically with NDSolve, this property is not respected. In particular

s[L_, c_] :=
NDSolve[{y''[x] == -Im[PolyLog[2, -E^(I*y[x])]], y'[
0] == c, y'[L] == c}, y, {x, 0, L}];
a = s[5, 0.2];
Evaluate[y[5] /. a] - Evaluate[y[0] /. a]


gives me an output $$5.02913$$, which is not nearly a multiple of $$2\pi$$. My guess is that something is wrong with the quadrature. Can you help me please? Maybe I have to use a special method, like a shooting method, or something? If so then what is the deal, are there any caveats?

Best, Kiryl

Replacing one IC seems to work:

ClearAll[f, F, inv, sol]
f[w_?NumericQ] := Im[PolyLog[2, -E^(I*w)]]
F[w_?NumericQ] := NIntegrate[f[t], {t, 0, w}]
inv[w_?NumericQ, wp_?NumericQ]:= 0.5*wp^2+F[w]

sol[l_, c_] := NDSolve[
{wp[x] == w'[x], wp'[x] == -f[w[x]], wp[ 0] == c, Mod[Abs[w[0]  - w[l]], 2*Pi] == 0},
{w, wp},
{x, 0, l},
MaxStepFraction -> 0.01,
MaxSteps -> Infinity,
Method -> {"FixedStep", Method -> Automatic}
];

out = sol[5, 0.2] // First ;
w[5] - w[0] /. out
{wp[0], wp[5]} /. out

F[w[5] /. out] - F[w[0] /. out]
Table[inv[w[x], wp[x]] /. out, {x, 0, 5, 0.5}]


I'd expect projection method to work here, but for some reason it throws an error:

sol[l_, c_] := NDSolve[
{wp[x] == w'[x], wp'[x] == -f[w[x]], wp[ 0] == c, wp[l] == c},
{w, wp},
{x, 0, l},
MaxStepFraction -> 0.01,
MaxSteps -> Infinity,
Method -> {"FixedStep", Method ->  {"Projection", Method -> Automatic,  "Invariants" ->{inv[w[x], wp[x]] }}}
];

out = sol[5, 0.2] // First ;
(* NDSolve::nnum1: The function value inv[w[0.],wp[0.]] is not a number when the arguments are {0.,{0.,0.}}.  *)


Edit

ClearAll[f, F, inv, sol]
f[w_?NumericQ] := Im[PolyLog[2, -E^(I*w)]]
F[w_?NumericQ] := NIntegrate[f[t], {t, 0, w}]
inv[w_?NumericQ, wp_?NumericQ] := 0.5*wp^2 + F[w]

ClearAll[solution] ;
Options[solution] = {MaxStepFraction -> 0.005, MaxSteps -> Infinity,
Method -> {"FixedStep",
Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 10}}} ;
solution[l_, c_,  opts : OptionsPattern[]] :=
NDSolve[{wp[x] == w'[x], wp'[x] == -f[w[x]], wp[0] == c,
Mod[Abs[w[0] - w[l]], 2*Pi] == 0}, {w, wp}, {x, 0, l}, opts] /;
c >= 0 ;
solution[l_, c_,  opts : OptionsPattern[]] :=
NDSolve[{wp[x] == w'[x], wp'[x] == -f[w[x]], wp[l] == c,
Mod[Abs[w[0] - w[l]], 2*Pi] == 0}, {w, wp}, {x, 0, l}, opts] /;
c < 0 ;

out = solution[5, 0.2] // First;
{wp[0], wp[5]} /. out
Table[inv[w[x], wp[x]] /. out, {x, 0, 5, 1}]
(* {0.2,0.19999537137167114} *)
(* \
{-1.6718455006855417,-1.6718458095149678,-1.671842921367624,-1.\
6718458808125778,-1.6718464107564595,-1.6718464264002912} *)

out = solution[5, -0.2] // First;
{wp[0], wp[5]} /. out
Table[inv[w[x], wp[x]] /. out, {x, 0, 5, 1}]
(* {-0.19999537137418164,-0.2000000000000036} *)
(* \
{-1.6718464264011181,-1.671846410757299,-1.6718458808139036,-1.\
6718429213688029,-1.671845809516107,-1.67184550068668} *)

• Thank you for the input, very insightful! The problem is that with this method, the derivative at $L$ is sometimes negative $c$ and it has to be equal at $x=0$ and $x=L$ Sep 6, 2023 at 23:03
• @KirylPesotski, indeed, see the edit, I still believe the project method to be a way to go, not sure why it is failing (guess it is related to the problem being bvp)
– I.M.
Sep 7, 2023 at 3:56