1
$\begingroup$

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

enter image description here

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}]

enter image description here

which is clearly not conserved.

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Bill Watts Dec 4 '18 at 1:49
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey Dec 4 '18 at 2:12
  • 2
    $\begingroup$ You also could try the Projection Method. $\endgroup$ – bbgodfrey Dec 4 '18 at 2:20
  • 5
    $\begingroup$ 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. $\endgroup$ – Chris K Dec 4 '18 at 3:17
4
$\begingroup$

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}]

fig1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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