Non Trivial Conserved Quantity Fails To Be Conserved After Solving System of ODE's Via NDSolve

Question

I'm trying to solve the following Darboux system of equations numerically.

As a result I have the following implementation of NDSolve

NDSolve[{w1'[t] == w2[t]*w3[t] - w1[t] (w2[t] + w3[t]), 
w2'[t] == w1[t]*w3[t] - w2[t] (w1[t] + w3[t]), 
w3'[t] == w1[t]*w2[t] - w3[t] (w1[t] + w2[t]), 
w1[0] == E^(2*1 - 2 - Sqrt[3]*3), w2[0] == E^(2*1 - 2 + Sqrt[3]*3), 
w3[0] == E^(2*1 + 2*2)}, {w1[t], w2[t], w3[t]}, {t, 0, 100}]

However when I solve this and plug it into one of four conserved quantities associated with the Darboux system I get

Plot[Re[N[
1/((Sqrt[2]
Sqrt[((-w1[t] + w2[t]) (w1[t] - w3[t]))/(w2[t] - 
w3[t])^2] ((-EllipticE[(w1[t] - w3[t])/(w2[t] - w3[t])] + 
EllipticK[(w1[t] - w3[t])/(w2[t] - w3[t])]) w2[t] + 
EllipticE[(w1[t] - w3[t])/(w2[t] - w3[t])] w3[
t]))/(\[Pi] Sqrt[((w1[t] - w2[t]) (w1[t] - w3[t]))/(
w2[t] - w3[t])])) /. s]], {t, 1, 100}]

which is clearly not conserved.

Does anyone know of a way to configure NDSolve so my solution leads to a perfectly conserved quantity? Thanks!

Answer 1

In numerical calculations, invariants are preserved with some accuracy. To improve accuracy, you can use a special method, for example

invariant = 
  1/((Sqrt[2] Sqrt[((-w1[t] + w2[t]) (w1[t] - w3[t]))/(w2[t] - 
            w3[t])^2] ((-EllipticE[(w1[t] - w3[t])/(w2[t] - w3[t])] + 
            EllipticK[(w1[t] - w3[t])/(w2[t] - w3[t])]) w2[t] + 
         EllipticE[(w1[t] - w3[t])/(w2[t] - w3[t])] w3[
           t]))/(\[Pi] Sqrt[((w1[t] - w2[t]) (w1[t] - w3[t]))/(w2[t] -
            w3[t])]));

s = NDSolve[{w1'[t] == w2[t]*w3[t] - w1[t] (w2[t] + w3[t]), 
    w2'[t] == w1[t]*w3[t] - w2[t] (w1[t] + w3[t]), 
    w3'[t] == w1[t]*w2[t] - w3[t] (w1[t] + w2[t]), 
    w1[0] == E^(2*1 - 2 - Sqrt[3]*3), 
    w2[0] == E^(2*1 - 2 + Sqrt[3]*3), w3[0] == E^(2*1 + 2*2)}, {w1[t],
     w2[t], w3[t]}, {t, 0, 100}, 
   Method -> {"Projection", Method -> "ExplicitRungeKutta", 
     "Invariants" -> invariant}, WorkingPrecision -> 50];


inv0 = Re[invariant /. First[s] /. t -> 0]

(*Out[]= 0.0753882728379682557368192924428895171708628882*)

Plot[Re[invariant /. First[s]] - inv0, {t, 1, 100}]