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

[Warren Weckesser]

On Tue, Oct 14, 2014 at 11:08 AM, Farid Afandiyev <faridafandiyev at gmail.com>
wrote:

> Hi Nicky,
>
> Correct, flow will stop after pressures are equalized. I would like it to
> stop integration at that point, i.e. after pressure equilibrium is reached
> odeint should stop calling derivative function.
>
>

Farid,

Your system of differential equations has an analytical solution.
To derive it, subtract the equation for P_2'(t) from the equation
for P_1'(t).  You'll end of with an equation of the form

   d(P_1(t) + P_2(t))/dt = -c * sqrt(P_1(t) - P_2(t))

(assuming P_1(0) > P_2(0), and where c is a constant expression of
the system parameters A, B, etc), or

   dx/dt = -c*sqrt(x)

where x(t) = P_1(t) - P_2(t).  That's a separable first order equation;
the solution is

  x(t) = (-c*t/2 + sqrt(x0))**2,

until x reaches 0 (i.e. until the vertex of the parabola).
x0 is the initial condition x(0) = P_1(0) - P_2(0).

If P_1(0) < P_2(0), you get

  x(t) = -(-c*t/2 + sqrt(-x0))**2

If you add Vc_1*B*P_1'(t) and Vc_2*B*P_2'(t), you find that

  Vc_1*B*P_1'(t) + Vc_2*B*P_2'(t) = 0

so

  Vc_1*P_1(t) + Vc_2*P_2(t) = constant = Vc_1*P_1(0) + Vc_2*P_2(0)

>From that and the solution x(t), you can get P_1(t) and P_2(t).

You can do the same type of calculation for the differential equations
for the two Va components (if you need them).

The attached script creates a plot comparing the analytical solution
to the solution found using odeint.  It generates the attached plot.


Warren




> Thanks,
>
> Farid.
>
> On Tue, Oct 14, 2014 at 5:05 PM, nicky van foreest <vanforeest at gmail.com>
> wrote:
>
>> Hi Farid,
>>
>> Please see below.
>>
>> [[  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]
>>>
>>
>> It appears that your P1 - P2 = 0 at this line. It might mean that the
>> pressure in both vessels is the same at that moment in time, and then the
>> flow should stop, right?  So, what should the integrator tell you after
>> this point in time?
>>
>> Nicky
>>
>>
>>>  [  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]]
>>>
>>>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20141015/4b73bbcd/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pressure_odeint_question.py
Type: text/x-python
Size: 3158 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20141015/4b73bbcd/attachment.py>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pressure.png
Type: image/png
Size: 73346 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20141015/4b73bbcd/attachment.png>