Scipy: Two ways of implementing a differential equation: two different solutions: answered

Question

I was trying to solve a differential equation for my thesis in chemistry and there I stumbled over a question regarding the differential equation solver "odeint" of scipy.

First I implemented the differential by the function CIDNP_1 (CIDNP is a chemical phenomena, that explains the unusual variables) according the the example on the scipy website. But the solution is way off even the right direction.

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate

R0 = 5e+5
kt = 5e5/R0
beta = 3/R0

def CIDNP_1(y, t):
    dP_dt, dQ_dt = y

    def R(t):
        return R0/(1 + kt*R0*t)

    dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2
    dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2
    return [dP_dt, dQ_dt]


def CIDNP_2(y, t):    
    dP_dt, dQ_dt = y

    def R(t):
        return R0/(1 + kt*R0*t)

    return [-kt*dP_dt*R(t) - kt*beta*(R(t))**2, \
            +kt*dP_dt*R(t) + kt*beta*(R(t))**2]

y0 = [-1, +1]
t = np.linspace(1e-9, 100e-6, 1e3)
sol_1 = scipy.integrate.odeint(CIDNP_1, y0, t)
sol_2 = scipy.integrate.odeint(CIDNP_2, y0, t)

Then I altered my solution to CIDNP_2 what gave the correct result but in my eyes there is no difference in the implementation as the variables dP_dt and dQ_dt are not altered in implementation CIDNP_1.

So as anyone can give me a hint why the implementation CIDNP_1 gives the wrong result I'd be very lucky as then at least the last two hours weren't completely lost.

Regards,

Jakob

Answer 1

In your first version you perform the updates not at the same time, as you execute the to lines

dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2
dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2

not simulatnously; therefore, you use the already updated dP_dt to update dQ_dt. This is a faulty implementation of the ODE system. Your second approach is better, as it does not have this kind of error. You either have to return the updated values directly, or you have to save the newly calculated values of dP_dt in another variable before calculating the new dQ_dt.

new_dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2
new_dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2

return [new_dP_dt, new_dQ_dt]

this would solve your problem.