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!

Improve this question
$\endgroup$
4
  • $\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$ Commented Dec 4, 2018 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$ Commented Dec 4, 2018 at 2:12
  • 2
    $\begingroup$ You also could try the Projection Method. $\endgroup$ Commented Dec 4, 2018 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$ Commented Dec 4, 2018 at 3:17

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks a lot Alex and to everyone who commented. I greatly appreciate it. $\endgroup$ Commented Dec 4, 2018 at 16:14

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.