Multiple Coupled Differential Equation solution in Python

Question

I have 4 ordinary differential equations that are coupled. The variables in the 4 equations are functions of time and space and one of them is second order in space.

\begin{equation} \frac{ \partial t_g}{ \partial t} = \frac {hp(t_m-tg)}{A(\rho_gCp_g)}- u\frac { \partial t_g}{ \partial x} \end{equation} \begin{equation} \frac{ \partial t_m}{ \partial t} = \frac{k_{sub}A_{sub}}{(\rho_mA_mCp_m)+(\rho_{sub}A_{sub}Cp_{sub})}\frac{ \partial^2 t_m}{ \partial x^2} + \frac{ h_mp(\rho_{vg}-\rho_{vm})h_{ads}}{ (\rho_mA_mCp_m)+(\rho_{sub}A_{sub}Cp_{sub})} - \frac {hp(t_m-t_g)}{(\rho_mA_mCp_m)+(\rho_{sub}A_{sub}Cp_{sub}} \end{equation} \begin{equation} \frac { \partial \rho_{vg}}{ \partial t} = - \frac {h_mp( \rho_{vg} -\rho_{vm})}{A} - \frac {u \partial \rho_{vg}}{ \partial x} \end{equation} \begin{equation} \frac { \partial \gamma_m}{ \partial t} = \frac {h_mp(\rho_{vg}-\rho_{vm})}{ \rho_mA_m} \end{equation} \begin{equation} \frac { \gamma_m}{ \gamma_{max}} = \frac {1}{1-c+ \frac{c}{ \phi}} \end{equation} \begin{equation} \phi = 4.09(10^{-9})t_m\rho_{vm}(e^\frac{-5196}{t_m}) \end{equation}

$t_g$, $t_m$, $\rho_{vg}$, $\rho_{vm}$, and $\gamma_m$ are the variables.

I am looking to solve this in Python. What is the best way to proceed?

I have the Boundary and Initial Conditions corresponding to these equations. But, having never used Python to solve such a problem, I am unsure whether to use a function like odeint or to type out the forward difference scheme.

Answer 1

You have a relation for $\gamma_m$: $$ \gamma_m = \frac {\alpha_7}{1-c+ \frac{c}{ \phi}} =: f(t_m, \rho_{vm}; t,x) $$ which can be used to get an equation of the form $$ \frac{\partial \gamma_m}{\partial t} = \beta_\rho(t_m, \rho_{vm})\,\frac{\partial \rho_{vm}}{\partial t} + \beta_t(t_m, \rho_{vm})\,\frac{\partial t_m}{\partial t} + \beta_c(t_m, \rho_{vm}) $$ You can then use the expression for the derivative of $\gamma_m$ $$ \frac { \partial \gamma_m}{ \partial t} = \alpha_6(\rho_{vg}-\rho_{vm}) $$ to write $$ \frac{\partial \rho_{vm}}{\partial t} = \alpha(t_m,\rho_{vg},\rho_{vm}) + \beta(t_m, \rho_{vm})\,\frac{\partial t_m}{\partial t} $$ Then the four equations that you have are $$ \begin{aligned} \frac{ \partial t_g}{ \partial t} &= \alpha_1\,(t_m-t_g)- u\frac { \partial t_g}{ \partial x} \\ \frac{ \partial t_m}{ \partial t} &= \alpha_2\frac{ \partial^2 t_m}{ \partial x^2} + \alpha_3(\rho_{vg}-\rho_{vm}) - \alpha_4(t_m-t_g)\\ \frac { \partial \rho_{vg}}{ \partial t} &= - \alpha_5( \rho_{vg} -\rho_{vm}) - u\frac {\partial \rho_{vg}}{ \partial x} \\ \frac{\partial \rho_{vm}}{\partial t} &= \alpha(t_m,\rho_{vg},\rho_{vm}) + \beta(t_m, \rho_{vm})\,\frac{\partial t_m}{\partial t} \end{aligned} $$ where $t_g$, $t_m$, $\rho_{vg}$, $\rho_{vm}$ are the variables you want to solve for.

Now imagine the spatial derivatives fixed so that these are ordinary differential equations with $t$ as the independent variable. You will still not be able to solve the system because only three of these equations are independent.

Answer 2

Just discretize the all the spatial term of every equation.

Here i'm concentrating only first equation other equation should be done in similar style. Lets use CD4 in space and RK4 in time because CD4 has good Dispersion preserving property and RK4 has relatively good stability. Please try to avoid Euler time integration at any cost if you are using central difference in space and the equation has dominant advection character because of stability issues.

Please note that we should use some lower order methods at boundary where we can't apply CD4 stencil else if the boundary is periodic we can apply CD4 over the whole domain.
\begin{equation} \frac{ \partial t_g}{ \partial t} = \frac {hp(({t_m}_i-{t_g}_i)}{A(\rho_gCp_g)}- u\frac { {t_g}_{i-2}-8*{t_g}_{i-1}+8{t_g}_{i+1}-{t_g}_{i+2}}{ 12*\delta x} \end{equation}

Form the initial condition you know $t_g$ and $t_m$ values in the domain at t=0 now we should march in time. Lets assume \begin{equation} F= \frac {hp(({t_m}_i-{t_g}_i)}{A(\rho_gCp_g)}- u\frac { {t_g}_{i-2}-8*{t_g}_{i-1}+8{t_g}_{i+1}-{t_g}_{i+2}}{ 12*\delta x} \end{equation}

This equation become

\begin{equation} \frac{ d t_g}{ \partial t} =F(t_g) \end{equation}

Using low storage RK4 method
this will become

\begin{equation} {t_g}^0 = {t_g}^n \end{equation} here you can substitute Initial condition value of $t_g$

${t_g}^1 = {t_g}^0 + \alpha_1 \delta t F({t_g}^0)$

${t_g}^2 = {t_g}^0+\alpha_2 \delta t F({t_g}^1)$

${t_g}^3 = {t_g}^0+\alpha_3 \delta t F({t_g}^2)$

${t_g}^4 = {t_g}^0 + \alpha_4 \delta t F({t_g}^3)$

${t_g}^{n+1}={t_g}^4 $

$\alpha_4=1/4$, $\alpha _3=1/3$, $\alpha _2=1/2$ and $\alpha _1=1$

Like this do this for other equation. Since we are not having any time derivative of $\rho_{vm}$ we should use those algebraic equations to calculate $\rho_{vm}$ using newly calculated other variables.

Please note that this kind of problem may cause dispersion and numerical diffusion so always use low CFL number -may be less than 0.5 and fine grid (this should be determined by trial and error if you know wave number of initial condition then it may be easy but most of the time we may not have single wave numbered initial condition).