# scipy odeint: sum of conservative ode equations does not remain zero as it is being solved — is this normal?

Let us say we have the following equations:

dy1/dt = f(y1, t)      [1]
dy2/dt = g(y2, t)      [2]


The equations are such that they are "conservative", i.e. the following condition should hold:

dy1/dt + dy2/dt = 0    [3]


Using scipy.odeint, I find that I can integrate conservative equations like this just fine for simple systems of ODEs.

However, for larger ones, I get the following issue.

Let's say that this is my derivative function:

def deriv_function(y0s, t):
...body defines equations 1, and 2...

print np.sum(ode)
return ode


Note the print statement.

I use scipy.odeint on deriv_function as follows:

odeint(deriv_fun, y0s, [0, 0.5])


Due to the print statement, the following gets printed out:

-1.38555833473e-13 <--- note, close to zero
-0.00679107743937
-0.0067907211796
-0.0135814423985
-0.0135810861584
-0.416522145214
-0.416523165887
-0.818209018574
-0.818211056221
-1.21864678558
-1.21864881584
-2.86735888212
-2.8673729885
-2.46855840934
-2.46856658088
-3.70632102566
-3.70631206163
-4.93200749506
-4.93200691488
-6.14577326158
-6.14577268283
-8.53799987128
-8.53799713959
-10.8839304356
-10.8839320212
-13.1845005689
-13.1845021725
-15.4406122011
-15.4406123927
-17.6531469917
-17.653147185
-24.6238415795
-24.6238498033
-31.1628867985
-31.1628947266
-37.2974784594
-37.2974547092
-35.463527103
-35.4635192949
-39.5777426955
-39.5777472677
-43.5137135424
-43.5137108017
-47.2791485087
-47.2791483993
-50.881424906
-50.8814244751
-54.3275507164
-54.3275502654    <--- note, not close to zero


For smaller systems of equations (not the same), the following gets printed out:

-1.13686837722e-13
0.0
0.0
0.0
-1.13686837722e-13
1.13686837722e-13
0.0
0.0
0.0
0.0
0.0
0.0
0.0
-1.13686837722e-13
0.0
0.0
0.0
1.13686837722e-13
0.0
0.0
0.0
0.0
-1.13686837722e-13
0.0
0.0
-1.13686837722e-13
0.0
0.0
0.0
0.0
0.0
0.0
-1.13686837722e-13
0.0
0.0
5.68434188608e-14
-5.68434188608e-14
-5.68434188608e-14
0.0


I have spent a lot of time trying to figure out if the issue has something to do with the way I have defined the equations, and I am now fairly certain that this is not the case. In order to confirm this, I would like to ask: is it possible at all for a non-conservative system of equations to first print an ode sum that is close to zero, but later print values that are not close to zero?

Another way to think of the question: in the larger system of equations, the ode sum initial prints out to be approximately zero, before increasing. What could be happening in the bowels of the solver for this problem to be occurring?

-

In general, numerical methods do not conserve quantities that may be conserved in a continuous system. The discretized system is typically non-conservative. A smaller time step (equivalently a smaller error tolerance) should help, but if you have a quantity that must be precisely conserved, you will have to use a method that is designed to conserve the desired quantity. Symplectic integrators are perhaps the most widely known and used examples of such a case.

As an example, consider a simple mass-spring oscillator:
$$\dfrac{dx}{dt} = v$$ $$\dfrac{dv}{dt} = -\dfrac{k}{m}x$$

The total energy $E = \frac{1}{2} \left( m v^2 + k x^2\right)$ is conserved since $$\dfrac{dE}{dt} = m v \dfrac{dv}{dt}+k x \dfrac{dx}{dt} = -kv x + kxv = 0$$

However, if you integrate this system using a Runge-Kutta type method the energy will not be conserved. A quick test using MATLAB's ode45 integrator gives the following energy plot which is not as dramatic as your case, but shows a clear decay in the energy:

This is a Hamiltonian system with Hamiltonian $H(x,v) = E(x,v) = \frac{1}{2} \left( m v^2 + k x^2\right)$ and would be a good candidate for the application of a symplectic integrator which would preserve the total energy.

-
Thanks for the explanation, and the pointer to what I should look into! I want to confirm though: what I am printing out is a symptom of the non-conservative effect of general numerical methods? –  user89 Jun 19 '14 at 4:59
As @DavidKetcheson said in his answer, it's impossible for me to rule out a bug in your code without a lot more information and testing. If it's not a bug, it's almost certainly due to the non-conservative nature of the numerical scheme you're using which is a form of discretization error. –  Doug Lipinski Jun 19 '14 at 12:24
You may want to try decreasing the tolerance or time step to see if the error decreases. If it does not, there is almost certainly a bug in you code. –  Doug Lipinski Jun 19 '14 at 12:25
@DougLipinski You said "symplectic integrators are perhaps the most widely known". What are other integrators that conserve particular quantities? –  Aurelius Jun 19 '14 at 13:09
@Aurelius My background is largely in fluid dynamics so I had in mind flux difference schemes that exactly conserve mass/momentum. There are many finite volume and finite difference schemes that meet this criteria. There are likely other examples as well, a search for conservative formulation turns up over 8 million results. –  Doug Lipinski Jun 19 '14 at 13:52

There are two possibilities:

1. You have a bug.
2. You are seeing an accumulation of discretization errors.

Without more detail, it is impossible for us to tell you which is actually happening in your simulation.

-