2
$\begingroup$

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

$\endgroup$
1

1 Answer 1

1
$\begingroup$

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`} *)
$\endgroup$
2
  • $\begingroup$ 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$ $\endgroup$ Sep 6, 2023 at 23:03
  • $\begingroup$ @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) $\endgroup$
    – I.M.
    Sep 7, 2023 at 3:56

Your Answer

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

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