Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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?

share|improve this question
add comment

2 Answers 2

up vote 5 down vote accepted

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:
energy decay

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.

share|improve this answer
    
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 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 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 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 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 at 13:52
show 4 more comments

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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