python求解微分方程

Question

我试图解决四个不同的微分方程。在谷歌搜索和研究后，我终于能够理解解决程序是如何工作的，但我无法具体地让这个问题正确运行。代码编译，但图形不正确。

我认为问题在于函数中的卷表达式，它将根据时间的变化而改变。这个体积在特定的时间将被用来求解微分方程的右手边。

时间向量的间隔、起点和终点是正确的。常量也是正确的。

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)

Answer 1

考虑到变量名称和方程结构的提示，您正在考虑一个化学反应。

A + B -> C + D

反应物a,b,c,d的浓度A,B,C,D有两个变化来源，

• 反应本身与反应速度，k*a*b和
• 反应物B在含浓度b0_in和体积速率V_in的溶液中的流入，导致V_in/V在所有组分中的相对浓度变化和V_in*b0_in/V在B中的加入。

这一切都很好地反映在你的系统的前四个方程式中。在处理卷的过程中，您正在以一种不一致的方式混合两种方法。或者V是t的一个已知函数，因此不是状态向量的一个组件，那么

V = V_reactor_initial + V_flow_inlet * t

或者将其视为状态的一个组成部分，则当前卷为

V = C[4]

体积变化率是

dVdt = V_flow_inlet.

修改第二种方法的代码如下所示

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])

有结果

​

​