Ask Your Question
0

Solving system of ODEs

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

spontarelliam gravatar image

updated 2017-08-01 05:09:27 -0600

FrédéricC gravatar image

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 ... (more)

edit retag flag offensive close merge delete

1 answer

Sort by » oldest newest most voted
0

answered 2016-04-25 16:00:09 -0600

updated 2016-04-25 20:53:45 -0600

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.

edit flag offensive delete link more

Your Answer

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

Add Answer

Question Tools

1 follower

Stats

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

Seen: 252 times

Last updated: Apr 25 '16