[SciPy-User] odeint - Excess work done on this call
Farid Afandiyev
faridafandiyev at gmail.com
Sat Oct 11 08:31:22 EDT 2014
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>
More information about the SciPy-User
mailing list