sympy dsolve returns incorrect answer

Question

I'm using sympy.dsolve to solve a simple ODE for a decay chain. The answer I get for different decay rates (e.g. lambda_1 > lambda_2) is wrong. After substituting C1=0, I get a simple exponential

-N_0*lambda_1*exp(-lambda_1*t)/(lambda_1 - lambda_2)

instead of the correct answer which has:

(exp(-lambda_1*t)-exp(-lambda_2*t)).

What am I doing wrong? Here is my code

sp.var('lambda_1,lambda_2 t')
sp.var('N_0')
daughter = sp.Function('N2',Positive=True)(t)

stage1 = N_0*sp.exp(-lambda_1*t)
eq = sp.Eq(daughter.diff(t),stage1*lambda_1 - daughter*lambda_2)

sp.dsolve(eq,daughter)

Answer 1

Your differential equation is (using shorter variable identifiers)

y' = A*N0*exp(-A*t) - B*y

Apply the integrating factor exp(B*t) to get the equivalent

(exp(B*t)*y(t))' = A*N0*exp((B-A)*t)

Integrate to get

exp(B*t)*y(t) = A*N0*exp((B-A)*t)/(B-A) + C

y(t) = A*N0*exp(-A*t)/(B-A) + C*exp(-B*t)

which is exactly what the solver computed.

Answer 2

Did you plug in C1 = 0, expecting to get a solution such that y(0) = 0? That's not how it works. C1 is an arbitrary constant in the formula, setting it to 0 is no guarantee that the expression will evaluate to 0 when t = 0.

Here is a correct approach, step by step

sol1 = sp.dsolve(eq, daughter)  

This returns a Piecewise because it's not known if two lambdas are equal:

Eq(N2(t), (C1 + N_0*lambda_1*Piecewise((t, Eq(lambda_2, lambda_1)), (-exp(lambda_2*t)/(lambda_1*exp(lambda_1*t) - lambda_2*exp(lambda_1*t)), True)))*exp(-lambda_2*t))

We can clarify that they are not:

sol2 = sol1.subs(Eq(lambda_2, lambda_1), False)

getting

Eq(N2(t), (C1 - N_0*lambda_1*exp(lambda_2*t)/(lambda_1*exp(lambda_1*t) - lambda_2*exp(lambda_1*t)))*exp(-lambda_2*t))

Next, we need C1 such that the right hand side turns into 0 when t = 0. So, take the right hand side, plug 0 for t, solve for C1:

C1 = Symbol('C1')
C1val = solve(sol2.rhs.subs(t, 0), C1, dict=True)[0][C1]

(It's not necessary to include dict=True but I like it, because it enforces uniform output of solve: it's an array of dictionaries.)

By the way, C1val is now N_0*lambda_1/(lambda_1 - lambda_2). Put it in:

sol3 = sol2.subs(C1, C1val).simplify()

and there you have it:

Eq(N2(t), N_0*lambda_1*(exp(lambda_1*t) - exp(lambda_2*t))*exp(-t*(lambda_1 + lambda_2))/(lambda_1 - lambda_2))

The expression is equivalent to N_0*lambda_1*(exp(-lambda_2*t) - exp(-lambda_1*t))/(lambda_1 - lambda_2) although SymPy seems loathe to combine the exponentials here.