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

Question

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?

Answer 1

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.

Answer 2

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.