Loading [MathJax]/jax/output/HTML-CSS/jax.js
Ask Your Question
0

Solving system of ODEs

asked 9 years ago

spontarelliam gravatar image

updated 7 years ago

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:

dwdt=QxVx(t)QyVηw(t)QzVηw(t) dxdt=QyVγy(t)+QzVζz(t)QxVx(t) dydt=QyVηw(t)QyVγy(t) dzdt=QzVηw(t)QzVζ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)

Preview: (hide)

1 Answer

Sort by » oldest newest most voted
0

answered 9 years ago

updated 9 years ago

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

[w(time)=L1(2(10g22272+407g2272+2409)161(g32272+89g22272+2511g2272+20286),g2272,time)+20161,x(time)=L1(2(31g22272+2276g2272+36303)161(g32272+89g22272+2511g2272+20286),g2272,time)+260161,y(time)=L1(6(2g22272+178g2272+1641)161(g32272+89g22272+2511g2272+20286),g2272,time)+12161,z(time)=L1(6(5g22272+445g2272+9657)161(g32272+89g22272+2511g2272+20286),g2272,time)+30161]

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.

Preview: (hide)
link

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: 9 years ago

Seen: 1,330 times

Last updated: Apr 26 '16