[SciPy-User] odeint - Excess work done on this call

[Farid Afandiyev]

Hello!

I'm newbie both in calculus and python/scipy so I apologize if this
question is too dumb.
I'm trying to model flow between two pressure vessels. Let's say we have
two points and a link between them like this
[Vc_1, P_1]=====A====[Vc_2, P_2]
Vc_1, Vc_2 are constants volumes of the nodes and P_1, P_2 varying pressure
of the nodes respectively.

I end up writing differential questions below. Never mind physical meaning
I just want to get math correct.

dP_1/dt = dVa/dt * 1/(Vc_1 * B)
dP_2/dt = dVa/dt * 1/(Vc_2 * B)

Here B is compressibility.

dVa/dt = A * K * sqrt(P_1 - P_2)

dVa is amount of "flow" between nodes.

K is some constant coefficient and A is link "weight". (P_1-P_2) can change
sign so I've adjusted this in the software.

Below is python program that I wrote to evaluate this.

#!/usr/bin/env python

import math
from scipy.integrate import odeint
from time import time
import numpy

B_compressibility = 0.0000033 #water compressibility
K = 0.747871759938 #coefficient

Vc_1 = 20
Vc_2 = 50
A = 0.01
P_1 = 4000
P_2 = 2000

def deriv(state, t):
    _P_1 = state[0]
    _P_2 = state[2]
    diff_P = _P_1 - _P_2
    flow_direction = math.copysign(1, diff_P)
    dVa = flow_direction * A * K * math.sqrt(abs(diff_P))
    dP_1 = -(dVa/Vc_1)/B_compressibility
    dP_2 =  (dVa/Vc_2)/B_compressibility
    return [dP_1, -dVa, dP_2, dVa]

if __name__ == '__main__':
    Va_1 = Vc_1 * P_1 * B_compressibility
    Va_2 = Vc_2 * P_2 * B_compressibility
    odeIterations = 10
    timeperiod = numpy.linspace(0.0,1.0, odeIterations)
    initial_state = [P_1, Va_1, P_2, Va_2]
    t0 = time()
    state_array = odeint(deriv, initial_state, timeperiod)
    t1 = time()
    print 'runtime %fs' %(t1-t0)
    state = state_array[odeIterations-1]
    print state_array
    P_1 = state[0]
    Va_1 = state[1]
    P_2 = state[2]
    Va_2 = state[3]

Below is output from program

Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.
runtime 0.157000s
[[  4.00000000e+03   2.64000000e-01   2.00000000e+03   3.30000000e-01]
 [  3.49242034e+03   2.30499743e-01   2.20303186e+03   3.63500257e-01]
 [  3.09580400e+03   2.04323064e-01   2.36167840e+03   3.89676936e-01]
 [  2.81015098e+03   1.85469965e-01   2.47593961e+03   4.08530035e-01]
 [  2.63546127e+03   1.73940444e-01   2.54581549e+03   4.20059556e-01]
 [  2.57173487e+03   1.69734501e-01   2.57130605e+03   4.24265499e-01]
 [  2.57142857e+03   1.69714286e-01   2.57142857e+03   4.24285714e-01]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]]
 lsoda--  at current t (=r1), mxstep (=i1) steps
       taken on this call before reaching tout
      In above message,  I1 =       500
      In above message,  R1 =  0.6110315150411E+00



Odeint gives correct results up to 7th line and then something goes
seriously wrong. I have searched in google and it looks like I'm not the
only one who struggles with scipy. Everybody suggests increasing mxstep but
that doesn't solve my problem. In addition it slows down the method
significantly. Somebody suggested reducing accuracy but I don't know how to
do that. Also I don't need super accuracy from the odeint. Couple of digits
after dot is more than enough for me. Any help is greatly appreciated!

Kind regards,

Farid.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20141011/70ab3721/attachment.html>