Skip to main content
Mathematica

Return to Answer

added 1464 characters in body
Source Link
I.M.
I.M.
  • 3.2k
  • 1
  • 14
  • 19

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`} *)

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`} *)
Source Link
I.M.
I.M.
  • 3.2k
  • 1
  • 14
  • 19

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.}}.  *)