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(vessels) and $P_1$, $P_2$ are varying pressures of the nodes respectively.
I end up writing differential questions below. Never mind physical meaning I just want to get math correct.
$\frac{\mathrm{dP_1} }{\mathrm{d} t} = \frac{\mathrm{dVa} }{\mathrm{d} t} \cdot \frac{1}{Vc_1 \cdot B}$
$\frac{\mathrm{dP_2} }{\mathrm{d} t} = \frac{\mathrm{dVa} }{\mathrm{d} t} \cdot \frac{1}{Vc_2 \cdot B}$
Here $B$ is compressibility.
$\frac{\mathrm{dVa} }{\mathrm{d} t} = A \cdot K \cdot \sqrt {P_1 - P_2}$
$\frac{\mathrm{dVa} }{\mathrm{d} t}$ is amount of "flow" or change of additional volume between nodes. $K$ is some constant coefficient and $A$ is link "throughput". ($P_1-P_2$) can change sign so I've adjusted for 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 #print 'IN ', state #print 'OUT ', [dP_1, -dVa, dP_2, dVa] 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) print state_array P_1, Va_1, P_2, Va_2 = state_array[odeIterations-1]
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.041000s [[ 4.00000000e+003 2.64000000e-001 2.00000000e+003 3.30000000e-001] [ 3.49242034e+003 2.30499743e-001 2.20303186e+003 3.63500257e-001] [ 3.09580400e+003 2.04323064e-001 2.36167840e+003 3.89676936e-001] [ 2.81015098e+003 1.85469965e-001 2.47593961e+003 4.08530035e-001] [ 2.63546127e+003 1.73940444e-001 2.54581549e+003 4.20059556e-001] [ 2.57173487e+003 1.69734501e-001 2.57130605e+003 4.24265499e-001] [ 2.57142857e+003 1.69714286e-001 2.57142857e+003 4.24285714e-001] [ 1.83357137e-299 1.80790662e-299 1.83145695e-299 1.87152166e-299] [ 1.83276935e-299 1.80681296e-299 1.83182150e-299 1.87141230e-299] [ 1.83379011e-299 1.80746916e-299 1.83320682e-299 1.83229543e-299]] 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. Decreasing accuracy is OK if that help as I don't need super accuracy from the odeint
. Couple of digits after dot is more than enough for me. Also I just need final values so in fact I want to decrease number of steps. Any help is greatly appreciated!