Chemical Reaction Differential Equations in Python

Concentrations on chemical species from mole balances are solved for 1, 2, and 4 species in Python with the Scipy.Integrate package ODEINT.

[@freekillerwhale]

Thank you so so so much, you saved my bachelor degree and sanity

[@Albert-PhotonChem]

I really appreciate the effort people like you put in sharing their knowledge. Very helpful tutorial. Thank you so much!!

[@apm]

I appreciate the feedback.

[@AJ-et3vf]

Awesome video sir! Thank you!

[@janjamniksrpcic4517]

Great tutorial!

[@ben5234]

Thanks this was really helpful

[@leonelamartes8829]

This is the best tutorial ever.

[@apm]

Wow! Thanks!

[@JoseAlvarez-dl3hm]

Hi great vids, could you handle in python this type of system of equations: diA/dt=La+diB/dt+iB(t) diB/dt=diA/dt+ Lb+iA(t) Because on all the videos I have seen about odeint in python I never saw a system which included a diferential in both sides of the equations. I truly appreciate any lead, thank you.

[@apm]

You are correct that Python ODEINT doesn't solve problems of this form. Here is a tutorial on typical problems: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations where you explicitly return the derivatives. Here are two options for your system of implicit differential equations: (1) use Python GEKKO that can handle implicit solutions such as you showed apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations or (2) perform a linear system solve to simultaneously solve for diA/dt and diB/dt within the ODEINT model function call. I recommend GEKKO because #2 may cause numerical problems with an implicit solve.

[@JoseAlvarez-dl3hm]

Thank you so much for the quick response, I will dive into that GEKKO option. Cheers!

[@oscargomez412]

Thank you very much for the video, very useful your channel. I have a system very similar to that of problem 3. In this case we have experimental values of the concentrations and we do not know the values of the kinetic constants "K's". We solved the ODEs assuming the values of the "K's". Any suggestion (video, exercise, etc) of how to adjust the values of these kinetic constants so that the result of the ODEs approaches the experimental results and thus determine the real "K's".

[@apm]

Here is a related problem with an estimated kinetic constant. I recommend the Gekko code with Example 3: apmonitor.com/do/index.php/Main/DynamicEstimation

[@oscargomez412]

@@apm Thank you very much for your prompt response, it is very helpful.

[@algerianinusa]

Thank you so much. Could you please include the noise, and solve it using Gillepsie Algorithm?

[@apm]

More information on dealing with noise is here: apmonitor.com/do/index.php/Main/DynamicData with estimators. Here is a similar example with noise added. apmonitor.com/pdc/index.php/Main/SimulateHIV The Gillepsie algorithm sounds interesting.

[@user-oh7on3mv2g]

Thanks a lot. But I got a question. I built two models to calculate total mass loss of reaction below, but why my second one was wrong? A -> kB + (1-k)gas def dmdT(m, T, A, E, n, k): m_plywood = m[0] m_residue = m[1] m_total = m[2] r1 = -A * exp(-E/R/T) * m_plywood**n / beta r2 = k * A * exp(-E/R/T) * m_plywood**n / beta dm_plywooddT = r1 dm_residuedT = r2 dm_totaldT = r1 + r2 return [dm_plywooddT, dm_residuedT, dm_totaldT] def dm_totaldT(m, T, A, E, n, k): r1 = -A * exp(-E/R/T) * m**n / beta r2 = k * A * exp(-E/R/T) * m**n / beta return r1 + r2

[@apm]

If m is size 3 then you need to return 3 derivatives. Here are some tutorials: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations with ODEINT and with Python Gekko: apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations

[@user-oh7on3mv2g]

@@apm Thanks a lot

[@rvilleg95]

Hello, I do the exact same code in Python and Jupyter and I get that the selectivity and the concentracion of C get to be negative, I have honestly checked the code line by line and haven't found the error. And there where you call the concentracion for the last exercise, when you say "conc=odeint(rxn,Z0,t)" wouldn't it be conc=odeint(rxn(Z0,t)) because you are passing data to a function? thanks again

[@apm]

+Rafael Villegas it is probably just an incorrect negative sign in your differential equation for the concentration of C (check equation 3 signs). Your notation makes sense but ODEINT requires them as 3 separate arguments. See docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html Great questions!

[@Antonio-lt1sp]

I tried to replicate this code in Anaconda / Spyder (Python 3.5), but I got a synthax error at "%matplotlib inline". What is this line of code doing, and why is it not recognized by Python 3.5-Spyder-Anaconda?

[@apm]

You can omit that line - you only need it for the IPython Notebook. You may need the plt.show() command at the end of your script to show the plot if you are using something other than the IPython Notebook.

[@Antonio-lt1sp]

Thank you!

[@user-kai516]

Hi, it's really a great video! As a beginner, I would like to ask: S=np.empty(len(cC)) for i in range(len(cC)): if abs(cC[i]+cD[i]>1e-10): S[i]=cC[i]/(cC[i]+cD[i]) else: S[i]=1.0 What does it mean?(problem 3) thanks a lot : )

[@apm]

The concentrations cC[0] and cD[0] are zero so this avoids divide by zero for the selectivity calculation.

[@user-kai516]

@@apm Okay, thanks for reply. I will keep learning from your videos 🙂

[@PiyushKumar-kf7ei]

How would you plot the conversion of that last reaction vs time?

[@apm]

Here is help on plotting: github.com/APMonitor/data_science/blob/master/04.%20Visualize.ipynb after you calculate conversion.

[@vengdasankumaresan9186]

Hi, it's really a great video! As a beginner, I would like to ask, y my ODEINT solver can't run but I follow ur steps last question. odeint() missing 2 required positional arguments: 'y0' and 't' is pops out. BTW the code is as follow: import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def rxn(Z, t): k1 = 1 k2 = 1.5 r1 = k1*Z[0]*Z[1] r2 = k2*Z[1]*Z[2] dAdt = -r1 dBdt = -r1 -r2 dCdt = r1 -r2 dDdt = r2 return[dAdt,dBdt,dCdt,dDdt] t = np.arange(0,3.01,0.2) Z0 = [1,1,0,0] Conc = odeint(rxn(Z0,t)) cA = Conc[:,0] cB = Conc[:,1] cC = Conc[:,2] cD = Conc[:,3] S = np.empty(len(cC)) for i in range(len(cC)): if abs(cC[i]+cD[i]>1e-10): S[i] = cC/(cC+cD) else: S[i] = 1.0 plt.plot(t,cA) plt.plot(t,cB) plt.plot(t,cC) plt.plot(t,cD) plt.xlabel(time) plt.ylabel(Concentration) plt.legend('Ca','Cb','Cc','Cd') I hope u can help me

[@apm]

Replace the odeint call with Conc = odeint(rxn,Z0,t) Here are additional examples that may help: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations

[@aditimehta4248]

Try this: import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def rxn(Z, t): k1 = 1 k2 = 1.5 r1 = k1*Z[0]*Z[1] r2 = k2*Z[1]*Z[2] dAdt = -r1 dBdt = -r1 -r2 dCdt = r1 -r2 dDdt = r2 return[dAdt,dBdt,dCdt,dDdt] t = np.arange(0,3.01,0.2) Z0 = [1,1,0,0] Z = odeint(rxn,Z0,t) cA = Z[:,0] cB = Z[:,1] cC = Z[:,2] cD = Z[:,3] S = np.empty(len(cC)) for i in range(len(cC)): if abs(cC[i]+cD[i]>1e-10): S[i] == cC/(cC+cD) else: S[i] = 1.0 plt.plot(t,cA,'r--') plt.plot(t,cB,'b-') plt.plot(t,cC,'g*') plt.plot(t,cD,'y-') plt.xlabel('time') plt.ylabel('Concentration') plt.legend(['cA','cB','cC','cD']) it works now.

[@vidyam89]

Can I do the same problem as discrete model? Please help!

[@apm]

You need to linearize and put into state space form as shown here apmonitor.com/pdc/index.php/Main/StateSpaceModel Once it is in state space form, you can convert to discrete linear form as shown here: docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.signal.StateSpace.html

[@devargyachakraborty9068]

Any suggestions how to solve a system of pde

[@apm]

Yes, here is some help. apmonitor.com/do/index.php/Main/PartialDifferentialEquations

[@NadineLynch10]

How can I get the worksheet for Jupyter?

[@apm]

See the download instructions at 0:25 or see problem 16 here: apmonitor.com/che263/index.php/Main/CourseHomework

[@NadineLynch10]

@@apm thanks I managed to figure it out late last night Lol but thanks!! What an awesome course/videos you have! It somehow sparked my coding instincts again, yesterday I couldn’t stop coding in python lol I even shared it with my friends

[@husseinalsajer4381]

thank you sooo fine please can you help me in this : x(2x+1)y '' +2(x+1)y' -2y =0 y'(1)= 0 y(3)-y'(3) = 31/9 y(x) = x+1+1/x this all question please if you can help me ?

[@apm]

Unfortunately, I can't help with specific questions but here are some resources in MATLAB: apmonitor.com/che263/index.php/Main/MatlabDynamicSim and Python: apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations