python solving differential equations

Question

I am attempting to solve four different differential equations. After googling and researching I was able to finally understand how the solver works but I can't get this problem specifically to run correctly. Code compiles but the graphs are incorrect.

I think the problem lies in the volume expression inside the function, which will change depending on how much time has passed. That volume at a specific time will then be used to solve the right hand side of the differential equations.

The intervals, starting point and ending point for the time vector is correct. Constants are also correct.

import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt

#defining all constants and initial conditions
k=2.2
CB0_inlet=.025
V_flow_inlet=.05
V_reactor_initial=5
CA0_reactor=.05
FB0=CB0_inlet*V_flow_inlet

def dcdt(C,t):
    #expression of how volume in reactor varies with time
    V=V_flow_inlet*t+C[4]   #C[4] is the initial reactor volume    ###we dont need things C to be C0 correct?

    #calculating right hand side of the four differential equations
    dadt=-k*C[0]*C[1]-((V_flow_inlet*C[0])/V)
    dbdt=((V_flow_inlet*(CB0_inlet-C[1]))/V)-k*C[0]*C[1]
    dcdt=k*C[0]*C[1]-((V_flow_inlet*C[2])/V)
    dddt=k*C[0]*C[1]-((V_flow_inlet*C[3])/V)

    return [dadt,dbdt,dcdt,dddt,V]  

#creating time array, initial conditions array, and calling odeint
t=np.linspace(0,500,100)
initial_conditions=[.05,0,0,0,V_reactor_initial]  # [CA0 CB0 CC0 CD0 
#V0_reactor]
C=integrate.odeint(dcdt,initial_conditions,t)
plt.plot(t,C)

I inlcuded them in my program just forgot to to do it here. sorry! — Casale
I thought your t vector was too sparse, but increasing the number of points from 100 to 10000 doesn't help. What seems most likely is that you've got a mistake somewhere in dcdt, as that's the most likely place where your integration could be going unstable, as signified by the exponential increase in the quantities you're calculating. — Perfi
Could you describe a little more on the context? V'(t)=V_in*t+V(t) indeed has an exponential solution V(t)=(V(0)+V_in)*exp(t)-(1+t)*V_in which might not be what you intended. — Lutz Lehmann

Lutz Lehmann Lutz Lehmann · Accepted Answer · 2018-04-07T12:11:35

Taking the hints of the variable names and equation structure, you are considering a chemical reaction

A + B -> C + D

There are 2 sources of changes in the concentration a,b,c,d of reactants A,B,C,D,

the reaction itself with reaction speed k*a*b and
the inflow of reactant B in a solution with concentration b0_in and volume rate V_in, which results in a relative concentration change of V_in/V in all components and an addition of V_in*b0_in/V in B.

This is all well reflected in the first 4 equations of your system. In the treatment of the volume, you are mixing two approaches in an inconsistent way. Either V is a known function of t and thus not a component of the state vector, then

V = V_reactor_initial + V_flow_inlet * t

or you treat it as a component of the state, then the current volume is

V = C[4]

and the rate of volume change is

dVdt = V_flow_inlet.

Modifying your code for the second approach looks like

import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt

#defining all constants and initial conditions
k=2.2
CB0_inlet=.025
V_flow_inlet=.05
V_reactor_initial=5
CA0_reactor=.05
FB0=CB0_inlet*V_flow_inlet

def dcdt(C,t):
    #expression of how volume in reactor varies with time
    a,b,c,d,V = C

    #calculating right hand side of the four differential equations
    dadt=-k*a*b-(V_flow_inlet/V)*a
    dbdt=-k*a*b+(V_flow_inlet/V)*(CB0_inlet-b)
    dcdt= k*a*b-(V_flow_inlet/V)*c
    dddt= k*a*b-(V_flow_inlet/V)*d

    return [dadt,dbdt,dcdt,dddt,V_flow_inlet]  

#creating time array, initial conditions array, and calling odeint
t=np.linspace(0,500,100)
initial_conditions=[.05,0,0,0,V_reactor_initial]  # [CA0 CB0 CC0 CD0 
#V0_reactor]
C=integrate.odeint(dcdt,initial_conditions,t)
plt.plot(t,C[:,0:4])

with the result

python solving differential equations

1 Answers