Deriving kinetic equations for a complex reaction system
In many cases, complex systems of even simple differential equations do not have simple and obvious solutions. Therefore, it could be too difficult and not reasonable to find fully analytical solutions.

Usually, numerical solutions can be easily calculated using iterative algorithms. However, such computer codes can be susceptible to calculation artefacts. Furthermore, the validation and verification of the model and the code can be a difficult issue.

I made a compromise and used an analytical approach when solving systems of differential equations, but I allowed the code to evaluate subexpressions involving constants numerically. Thus, I made numerical calculations to combine all constant parameter values as soon as possible in order to limit the number of parameters in the equations (using excessive digits in order to avoid accumulation of round-off errors).

I found that this approach can be easily and quickly implemented using the commercial software Maple. However, other similar mathematical software could be suitable, too.

I successfully used this approach for the evaluation of radioactive decay of the inventory of a nuclear reactor core in dependence of time after shutdown. The calculation of the radioactive decay of the relevant radionuclides involved more than 100 individual differential equations, which mainly describe chains of simple first-order reactions, e.g.:

$$\frac{\mathrm da_{\text{I-131}}(t)}{\mathrm dt}=-\lambda_{\text{I-131}}\cdot a_{\text{I-131}}(t)+0.778\cdot\lambda_{\text{I-131}}\cdot a_{\text{Te-131m}}(t)+1\cdot\lambda_{\text{I-131}}\cdot a_{\text{Te-131}}(t)$$

&gt; de_a_I131:=diff(a_I131(t),t)=-l_I131*a_I131(t)+0.778*l_I131*a_Te131m(t)+1*l_I131*a_Te131(t);


A large number of further differential equations could be added in order to model transport phenomena during postulated accidents (e.g. loss of coolant).

Differential equations that were dependent on each other and the corresponding initial conditions were grouped into systems of differential equations, e.g.:

&gt; sys_131:={de_a_Xe131m,de_a_I131,de_a_Te131,de_a_Te131m,a_Xe131m(t0)=core_Xe131m,a_I131(t0)=core_I131,a_Te131(t0)=core_Te131,a_Te131m(t0)=core_Te131m};


The systems of differential equations were solved analytically and the subexpressions involving constant parameter values were evaluated numerically.

&gt; sol_131:=evalf(value(dsolve(sys_131,[a_Xe131m(t),a_I131(t),a_Te131(t),a_Te131m(t)])));


The results could be easily rearranged into a suitable form, e.g.

$$\begin{aligned}a_\text{I-131}(t) ={} &amp; -1.6030637832147524320\cdot10^{17}\cdot\mathrm e^{-6.4180294496291232354\cdot10^{-6}\cdot t}\\&amp;-7.9769363886868464215\cdot10^{15}\cdot\mathrm e^{-4.6209812037329687295\cdot10^{-4}\cdot t}\\&amp;+4.7213833147101620896\cdot10^{18}\cdot\mathrm e^{-9.9782796169607015473\cdot10^{-7}\cdot t}\end{aligned}$$

The results could also be directly used for further calculations within the same worksheet.