python - How to solve a system of differential equations using scipy.odeint

Question

I want to solve a system of equations using odeint and I get the following error:

  File "C:", line 45, in <module>
    C_B = odeint(dC_Bdt,C_B0,t)

  File "C:Anaconda3envsChemEnglibsite-packagesscipyintegrateodepack.py", line 233, in odeint
    int(bool(tfirst)))

RuntimeError: The array return by func must be one-dimensional, but got ndim=2.

My code is:

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

# Create time domain
t = np.linspace(0,100,100)

# Parameters 
A = 1*10**(13) # Arrhenius constant
T = 293.15 # Temperature [K]
E_a= 80000 # Activation energy [J/mol]
R = 8.31 # Ideal gas constant [J/molK]
rho = 1000 # density [kg/m3]
F_in = 0.2 # Inlet flowrate [m3/s]
h = 2.1 # Height of reactor
A_= 1 # Cross-sectional area of reactor [m2]


# Initial condition
C_A0 = 3 # Initial concentration [mol/m3]
C_B0 = 0 # Initial concentration [mol/m3]
m0 = 0 # Initial mass in tank [kg]

def dmdt(F_out,t): # Mass balance
    return rho*(F_in-F_out)

def dC_Adt(C_A,t): # Concentration balance for A
    dC_Adt = (F_in*C_A-F_out*C_A)/V-k*C_A
    return dC_Adt

def dC_Bdt(C_B,t): # Concentration balance for B
    dC_Bdt = (F_in*C_B-F_out*C_B)/V+k*C_A
    return dC_Bdt #<-- needs to be 1D but is 2D


V = A_*h # Reactor volume [m3]
k=A*np.exp(-E_a/(R*T)) # Reaction rate constant
F_out = F_in # Steady state
dmdt = odeint(dmdt,m0,t)
C_A = odeint(dC_Adt,C_A0,t)
C_B = odeint(dC_Bdt,C_B0,t)

# Plot
plt.figure()
plt.plot(t,C_A,'b-',label='C_A')
plt.plot(t,C_B,'r--',label='C_B')
plt.legend(loc='best')
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Concentration [mol/m3]')

A similar question was asked at: How to fix the error: The array return by func must be one-dimensional, but got ndim=2

But the answer to that was to change the order of the function's dependencies. That doesn't work for me.

I think the problem arises from the fact that I use an output of a different function (C_A), to evaluate C_B. So could it be that dC_Bdt somehow has C_A and C_B as ouputs or at least links them somehow which is why when using odeint I get that error?

Not sure what to do from here? :( Many Thanks

Answer 1

I modified your code, in order to solve the three coupled equations:

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

# Parameters 
A = 1e13 # Arrhenius constant
T = 293.15 # Temperature [K]
E_a = 80000 # Activation energy [J/mol]
R = 8.31 # Ideal gas constant [J/molK]
rho = 1000 # density [kg/m3]
F_in = 0.2 # Inlet flowrate [m3/s]
h = 2.1 # Height of reactor
A_ = 1 # Cross-sectional area of reactor [m2]

V = A_*h  # Reactor volume [m3]
k = A*np.exp(-E_a/(R*T)) # Reaction rate constant
F_out = F_in # Steady state

def dUdt(U, t):
    m, C_A, C_B = U

    dmdt = rho*(F_in - F_out)
    dC_Adt = ( F_in*C_A - F_out*C_A )/ V-k*C_A
    dC_Bdt = ( F_in*C_B - F_out*C_B )/ V+k*C_A

    return [dmdt, dC_Adt, dC_Bdt]

# Create time domain
t_span = np.linspace(0, 100, 30)

# Initial condition
C_A0 = 3 # Initial concentration [mol/m3]
C_B0 = 0 # Initial concentration [mol/m3]
m0 = 0 # Initial mass in tank [kg]

Uzero = [m0, C_A0, C_B0]

solution = odeint(dUdt, Uzero, t_span)

# plot
plt.plot(t_span, solution[:, 0], label='masse');
plt.plot(t_span, solution[:, 1], label='C_A');
plt.plot(t_span, solution[:, 2], label='C_B');
plt.legend();
plt.xlabel('time');