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

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 == E^(2*1 - 2 - Sqrt*3), w2 == E^(2*1 - 2 + Sqrt*3),
w3 == 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
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!

• Don't know about your conserved quantities, but your code gives results that satisfy each differential equation and all initial conditions which you can see by plugging in values for t. Dec 4, 2018 at 1:49
• Numerical computations often have roundoff errors, and this is the case here. If you wish the relative error to be smaller use a larger WorkingPrecision. Dec 4, 2018 at 2:12
• You also could try the Projection Method. Dec 4, 2018 at 2:20
• Perhaps this is obvious, but notice that the plot's y-axis ranges from 0.0753880 to 0.075883, so the quantity is almost conserved. Dec 4, 2018 at 3:17

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

invariant =
1/((Sqrt 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 == E^(2*1 - 2 - Sqrt*3),
w2 == E^(2*1 - 2 + Sqrt*3), w3 == 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}] • Thanks a lot Alex and to everyone who commented. I greatly appreciate it. Dec 4, 2018 at 16:14