# Multiple Coupled Differential Equation solution in Python

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.

$$\frac{ \partial t_g}{ \partial t} = \frac {hp(t_m-tg)}{A(\rho_gCp_g)}- u\frac { \partial t_g}{ \partial x}$$ $$\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}}$$ $$\frac { \partial \rho_{vg}}{ \partial t} = - \frac {h_mp( \rho_{vg} -\rho_{vm})}{A} - \frac {u \partial \rho_{vg}}{ \partial x}$$ $$\frac { \partial \gamma_m}{ \partial t} = \frac {h_mp(\rho_{vg}-\rho_{vm})}{ \rho_mA_m}$$ $$\frac { \gamma_m}{ \gamma_{max}} = \frac {1}{1-c+ \frac{c}{ \phi}}$$ $$\phi = 4.09(10^{-9})t_m\rho_{vm}(e^\frac{-5196}{t_m})$$

$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.

• Welcome to SciComp.SE. Are you sure that you have a system of ODEs and not PDEs? Can you add your equations to the question? Dec 3, 2015 at 0:08
• (∂t_g)/∂t=(hp*(t_m- t_g))/(Aρ_g* Cp_g ) - (u∂t_g)/∂x Dec 3, 2015 at 15:06
• (∂t_m)/∂t = (k_subA_sub)/((ρ_m A_mCp_m) + (ρ_sub A_subCp_sub))*(∂^2 t_m)/(∂x^2 ) + (h_mpρ_vg- ρ_vm)*h_ads)/((ρ_m A_mCp_m)+(ρ_sub A_subCp_sub ) - (hp*(t_m-) t_g))/((ρ_m* A_mCp_m)+(ρ_sub A_sub*Cp_sub)) Dec 3, 2015 at 15:10
• Please, add them to the question. You can use LaTeX syntax to format them. Dec 3, 2015 at 15:30
• You have 5 unknowns and 4 equations. The system is not an ODE as you stated before... Dec 4, 2015 at 0:50

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.

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.
$$\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}$$

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 $$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}$$

This equation become

$$\frac{ d t_g}{ \partial t} =F(t_g)$$

Using low storage RK4 method
this will become

$${t_g}^0 = {t_g}^n$$ 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).

• Sorry right now I'm not having any sample code. I shall improve this answer when i get free time
– AGN
Dec 8, 2015 at 5:57
• Thanks for the tip on discretization !! That really helped . Won't the equation be a function of $t_g$ and $t_m$ ? And if so , both will vary with space and time again . The Initial condition is at t=0, at x=0 . To do a CD4 method discretization will the initial condition we have now be sufficient ? Dec 10, 2015 at 17:27
• Sorry i forget to mention, u could use lower order near boundary,where CD4 stencil can't be applied. For the first equation initial condition supposed to be $t_g(x,0)=f(x)$, Second equation requires two B.C and one initial condition. If initial condition is provided like $t_g(0,0) =c$. We should assume some condition except boundary. In that case we are interested in only steady state.
– AGN
Dec 11, 2015 at 1:27
• @Sanjeev For the first equation $t_m$ acts like source term. $t_g$ acts like source term in second equation. that I have discretized in $F$
– AGN
Dec 11, 2015 at 1:30