# 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[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!

• 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. – Bill Watts Dec 4 '18 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. – bbgodfrey Dec 4 '18 at 2:12
• You also could try the Projection Method. – bbgodfrey Dec 4 '18 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. – Chris K Dec 4 '18 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[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}]


• Thanks a lot Alex and to everyone who commented. I greatly appreciate it. – Daniel Berkowitz Dec 4 '18 at 16:14