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(time)=Vγm0ζ(ηγ+(ηγ+η)ζ)V+γζ+L−1(−V4g22562γm0ζ+(QxQz−Q2z)V2γ2m0ζ2−(QxQzηm0ζ2+(Q2x−QxQz)ηγ2m0)V3−((Qxηγm0+Qxηm0ζ)V4−((Qx−Qz)γ2m0ζ+Qzγm0ζ2)V3)g2562(V3g32562+(QxV3η+((Qx−Qz)γ+Qzζ+Qx)V2)g22562+(Q2xQz−QxQ2z)γζ+((Q2xQz−QxQ2z)ηγ+((Q2xQz−QxQ2z)ηγ+(Q2xQz−QxQ2z)η)ζ)V+((QxQzηζ+Q2xη+(Q2x−QxQz)ηγ)V2+((Q2x−QxQz)γ+(QxQz+(QxQz−Q2z)γ)ζ)V)g2562)((ηγ+(ηγ+η)ζ)V+γζ),g2562,time) x(time)=Vηγm0ζ(ηγ+(ηγ+η)ζ)V+γζ+L−1((QxQz−Q2z)Vγ2m0ζ2+(QxQzη2m0ζ2+(Q2x−QxQz)η2γ2m0)V3+((QxQz−Q2z)ηγ2m0ζ+(QxQz−Q2z)ηγm0ζ2)V2+(V3γm0ζ+(ηγm0+ηm0ζ)V4)g22562+((Qxη2γm0+Qxη2m0ζ)V4+((Qx−Qz)ηγ2m0+Qxηγm0ζ+Qzηm0ζ2)V3+((Qx−Qz)γ2m0ζ+Qzγm0ζ2)V2)g2562(V3g32562+(QxV3η+((Qx−Qz)γ+Qzζ+Qx)V2)g22562+(Q2xQz−QxQ2z)γζ+((Q2xQz−QxQ2z)ηγ+((Q2xQz−QxQ2z)ηγ+(Q2xQz−QxQ2z)η)ζ)V+((QxQzηζ+Q2xη+(Q2x−QxQz)ηγ)V2+((Q2x−QxQz)γ+(QxQz+(QxQz−Q2z)γ)ζ)V)g2562)((ηγ+(ηγ+η)ζ)V+γζ),g2562,time) y(time)=Vηm0ζ(ηγ+(ηγ+η)ζ)V+γζ+L−1(−V4ηg22562m0ζ+(QxQzη+(QxQz−Q2z)ηγ)V2m0ζ2+(QxQzη2m0ζ2+QxQzη2m0ζ−(Q2x−QxQz)η2γm0)V3+(QxV4η2m0ζ+(Qzηm0ζ2+((Qx−Qz)ηγ+Qxη)m0ζ)V3)g2562(V3g32562+(QxV3η+((Qx−Qz)γ+Qzζ+Qx)V2)g22562+(Q2xQz−QxQ2z)γζ+((Q2xQz−QxQ2z)ηγ+((Q2xQz−QxQ2z)ηγ+(Q2xQz−QxQ2z)η)ζ)V+((QxQzηζ+Q2xη+(Q2x−QxQz)ηγ)V2+((Q2x−QxQz)γ+(QxQz+(QxQz−Q2z)γ)ζ)V)g2562)((ηγ+(ηγ+η)ζ)V+γζ),g2562,time) z(time)=Vηγm0(ηγ+(ηγ+η)ζ)V+γζ+L−1(−V4ηg22562γm0−(QxQzη2m0ζ−((Q2x−QxQz)η2γ2+(Q2x−QxQz)η2γ)m0)V3+((QxQz−Q2z)ηγ2m0ζ+(Q2x−QxQz)ηγ2m0)V2+(QxV4η2γm0+(Qzηγm0ζ+((Qx−Qz)ηγ2+Qxηγ)m0)V3)g2562(V3g32562+(QxV3η+((Qx−Qz)γ+Qzζ+Qx)V2)g22562+(Q2xQz−QxQ2z)γζ+((Q2xQz−QxQ2z)ηγ+((Q2xQz−QxQ2z)ηγ+(Q2xQz−QxQ2z)η)ζ)V+((QxQzηζ+Q2xη+(Q2x−QxQz)ηγ)V2+((Q2x−QxQz)γ+(QxQz+(QxQz−Q2z)γ)ζ)V)g2562)((ηγ+(ηγ+η)ζ)V+γζ),g2562,time)
Can anyone explain what I've done wrong?