Ask Your Question

# Solving system of ODEs

I'm trying to solve a system of 4 ODEs using Sage. I'm able to get a solution, but it's not at all what I expect. Here are the 4 equations:

$$\frac{dw}{dt} = \frac{Q_x}{V} x(t) - \frac{Q_y}{V}\eta w(t) - \frac{Q_z}{V} \eta w(t)$$ $$\frac{dx}{dt} = \frac{Q_y}{V} \gamma y(t) + \frac{Q_z}{V} \zeta z(t) - \frac{Q_x}{V}x(t)$$ $$\frac{dy}{dt} = \frac{Q_y}{V} \eta w(t) - \frac{Q_y}{V} \gamma y(t)$$ $$\frac{dz}{dt} = \frac{Q_z}{V} \eta w(t) - \frac{Q_z}{V} \zeta z(t)$$

Here is how I've coded the problem in Sage.

eta, zeta, gamma = var('eta, zeta, gamma')
m_0, Qx, Qy, Qz, V = var('m_0, Q_x, Q_y, Q_z, V')
t = var('time')
w=function('w', t)
x=function('x', t)
y=function('y', t)
z=function('z', t)
dw = diff(w,t) - Qx/V*x + Qy/V*eta*w + Qz/V*eta*w == 0
dx = diff(x,t) - Qy/V*gamma*y - Qz/V*zeta*z + Qx/V*x == 0
dy = diff(y,t) - Qy/V*eta*w + Qy/V*gamma*y == 0
dz = diff(z,t) - Qz/V*eta*w + Qz/V*zeta*z == 0
eqs = [dw, dx, dy, dz]
vars = [w, x, y, z]
ics = [0, 0, m_0, 0, 0]
ans = desolve_system(eqs, vars, ics, ivar=t)
print ans


Where m_0 is an initial mass that exists in the x volume at time 0. All other volumes are empty at time 0.

I'm expecting relatively simple exponential equations as the solution, but what I get in return is some significantly more complex.

$$w\left(\mathit{time}\right) = \frac{V \gamma m_{0} \zeta}{{\left(\eta \gamma + {\left(\eta \gamma + \eta\right)} \zeta\right)} V + \gamma \zeta} + \mathcal{L}^{-1}\left(-\frac{V^{4} g_{2562}^{2} \gamma m_{0} \zeta + {\left(Q_{x} Q_{z} - Q_{z}^{2}\right)} V^{2} \gamma^{2} m_{0} \zeta^{2} - {\left(Q_{x} Q_{z} \eta m_{0} \zeta^{2} + {\left(Q_{x}^{2} - Q_{x} Q_{z}\right)} \eta \gamma^{2} m_{0}\right)} V^{3} - {\left({\left(Q_{x} \eta \gamma m_{0} + Q_{x} \eta m_{0} \zeta\right)} V^{4} - {\left({\left(Q_{x} - Q_{z}\right)} \gamma^{2} m_{0} \zeta + Q_{z} \gamma m_{0} \zeta^{2}\right)} V^{3}\right)} g_{2562}}{{\left(V^{3} g_{2562}^{3} + {\left(Q_{x} V^{3} \eta + {\left({\left(Q_{x} - Q_{z}\right)} \gamma + Q_{z} \zeta + Q_{x}\right)} V^{2}\right)} g_{2562}^{2} + {\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \gamma \zeta + {\left({\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta \gamma + {\left({\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta \gamma + {\left(Q_{x}^{2} Q_{z} - Q_{x} Q_{z}^{2}\right)} \eta\right)} \zeta\right)} V + {\left({\left(Q_{x} Q_{z ... edit retag close merge delete ## 1 answer Sort by » oldest newest most voted The answer returned consists of initial constants plus inverse Laplace transforms of fractions with quadratic numerators and cubic denominators. This is easier to see if you assign numeric values to all coefficients. For example eta=3; zeta=2; gamma=5 m_0=2; Qx=3; Qy=7; Qz=6; V=1 t = var('time') w=function('w', t) x=function('x', t) y=function('y', t) z=function('z', t) dw = diff(w,t) - Qx/V*x + Qy/V*eta*w + Qz/V*eta*w == 0 dx = diff(x,t) - Qy/V*gamma*y - Qz/V*zeta*z + Qx/V*x == 0 dy = diff(y,t) - Qy/V*eta*w + Qy/V*gamma*y == 0 dz = diff(z,t) - Qz/V*eta*w + Qz/V*zeta*z == 0 eqs = [dw, dx, dy, dz] vars = [w, x, y, z] ics = [0, 0, m_0, 0, 0] ans = desolve_system(eqs, vars, ics, ivar=t) show( ans )  produces the answer$$ \left[w\left(\mathit{time}\right) = \mathcal{L}^{-1}\left(-\frac{2 {\left(10 g_{2272}^{2} + 407 g_{2272} + 2409\right)}}{161 {\left(g_{2272}^{3} + 89 g_{2272}^{2} + 2511 g_{2272} + 20286\right)}}, g_{2272}, \mathit{time}\right) + \frac{20}{161}, x\left(\mathit{time}\right) = \mathcal{L}^{-1}\left(\frac{2 {\left(31 g_{2272}^{2} + 2276 g_{2272} + 36303\right)}}{161 {\left(g_{2272}^{3} + 89 g_{2272}^{2} + 2511 g_{2272} + 20286\right)}}, g_{2272}, \mathit{time}\right) + \frac{260}{161}, y\left(\mathit{time}\right) = \mathcal{L}^{-1}\left(-\frac{6 {\left(2 g_{2272}^{2} + 178 g_{2272} + 1641\right)}}{161 {\left(g_{2272}^{3} + 89 g_{2272}^{2} + 2511 g_{2272} + 20286\right)}}, g_{2272}, \mathit{time}\right) + \frac{12}{161}, z\left(\mathit{time}\right) = \mathcal{L}^{-1}\left(-\frac{6 {\left(5 g_{2272}^{2} + 445 g_{2272} + 9657\right)}}{161 {\left(g_{2272}^{3} + 89 g_{2272}^{2} + 2511 g_{2272} + 20286\right)}}, g_{2272}, \mathit{time}\right) + \frac{30}{161}\right] 

A fraction with a cubic denominator has a partial fraction decomposition that in general is three fractions with constant numerators and linear denominators. The constants that appear in the decomposition may be real or complex, depending on the values of the input coefficients. The inverse Laplace transform of each of these three fractions is an exponential or a circular function, depending on the values of the constants.

Your answer does consist of exponentials. It looks complicated because the 4x4 matrix of coefficients for the system, while relatively simple, leads to combinations of coefficients that are not simple, so that the inverse Laplace transforms cannot be evaluated in general form.

more

## Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

## Stats

Asked: 2015-07-02 13:52:10 -0600

Seen: 252 times

Last updated: Apr 25 '16