Master Scientific Visualization Using ParaView

Radioactive Decay Chain Simulator Using Bateman Equations and Monte Carlo Methods in Matlab

Author : Waqas Javaid

Abstract

This article presents a comprehensive and numerically stable simulator for modeling radioactive decay chains using analytical, numerical, and stochastic approaches. The framework combines the Bateman analytical solution with stiff and non-stiff ordinary differential equation solvers to accurately track nuclide populations over time. A transition matrix formulation is used to ensure physically consistent decay behavior and guaranteed non-negative solutions [1]. In addition to deterministic modeling, a Monte Carlo decay engine is implemented to capture probabilistic decay events and validate analytical predictions [2]. The simulator computes nuclide activities, total decay rates, and energy release profiles, enabling detailed physical interpretation. Sensitivity analysis is included to quantify how variations in decay constants affect final isotope populations. Multiple visualization outputs such as population curves, activity plots, solver comparisons, and sensitivity maps support deeper insight into decay dynamics [3]. The method is robust for both stiff and long half-life chains and avoids common numerical instabilities. This integrated radioactive decay chain simulation framework is suitable for research, teaching, and nuclear system analysis [4]. It provides an accurate and extensible foundation for advanced nuclear decay modeling and computational physics studies.

  1. Introduction

Radioactive decay chains play a fundamental role in nuclear physics, radiation engineering, geochronology, and medical isotope analysis, where unstable nuclides transform into daughter products through sequential decay processes.

Figure 1: Radioactive Decay Chain Simulator  Analytical, Numerical, and Monte Carlo Modeling of Nuclear Transmutation and Energy Release

Understanding how these chains evolve over time is essential for predicting isotope populations, radiation activity, and total energy release in nuclear systems. Because decay chains often include multiple isotopes with widely varying half-lives, their mathematical description leads to coupled differential equations that can be stiff and numerically challenging [5]. The classical Bateman equations provide an analytical solution for linear decay chains and remain a cornerstone for theoretical analysis, but practical implementations must address numerical stability and floating-point limitations [6]. Modern computational approaches complement analytical formulas with numerical ODE solvers that can efficiently handle stiff and non-stiff regimes. In addition, stochastic Monte Carlo simulation offers a particle-level interpretation by modeling individual decay events probabilistically, providing validation and physical intuition [7]. A robust decay chain simulator should therefore integrate analytical, deterministic, and stochastic techniques within a single framework. Such integration enables cross-verification of results and improves confidence in predictions [8].

Table 1: Final Activities

NuclideActivity (decays/s)
N11.159e-47
N23.694e-11
N33.544e-02
N46.037e-01
Stable0.000e+00

Accurate computation of activity and energy release further supports applications in reactor analysis, shielding design, and radiation safety assessment. Sensitivity analysis is also important because small uncertainties in decay constants can significantly affect long-term isotope distributions. Visualization of population evolution, activity curves, and solver behavior enhances interpretability for both researchers and students. Ensuring non-negative populations and physically consistent transitions is critical for credible simulations. Numerically stable formulations help avoid overflow, cancellation, and instability errors in long time horizons. Computational tools now allow high-resolution time grids and efficient matrix formulations for decay operators [9]. These advances make it feasible to build comprehensive and reliable decay chain models. This work presents a unified radioactive decay chain simulation framework that combines Bateman analytical solutions, stiff and non-stiff ODE solvers, and Monte Carlo methods [10]. The approach emphasizes numerical safety, physical correctness, and extensibility. It supports activity tracking, energy release estimation, and parameter sensitivity mapping. The resulting simulator serves as a practical platform for nuclear science education and applied computational physics.

1.1 Background and Importance of Radioactive Decay Chains

Radioactive decay chains are fundamental processes in nuclear science where unstable parent nuclides transform into a sequence of daughter isotopes until a stable product is formed. These chains appear in natural radioactive series, nuclear reactors, medical isotope production, and environmental radiation studies. Accurate modeling of decay chains is essential for predicting isotope concentrations and radiation levels over time [11]. The time evolution of these systems directly affects safety assessments, shielding design, and waste management strategies. In geoscience, decay chains are used for age dating of rocks and archaeological materials. In medicine, controlled decay chains support diagnostic imaging and targeted radiotherapy. Because each nuclide has its own half-life and decay constant, the overall system exhibits multi-scale time behavior [12]. This creates both analytical and computational challenges. Reliable simulation tools help researchers and engineers understand these coupled transformations. Therefore, decay chain modeling remains an important topic in computational nuclear physics.

1.2 Mathematical Modeling Foundations

The mathematical description of a radioactive decay chain is based on coupled first-order differential equations derived from decay laws. Each equation represents the rate of change of a nuclide population as a balance between its own decay and production from its parent. The Bateman equations provide a classical analytical solution for linear decay chains with known decay constants [13]. These closed-form expressions are elegant but can become numerically unstable when decay constants are close in value. Matrix formulations offer another perspective by representing the chain as a linear dynamical system. This approach enables compact representation and compatibility with numerical solvers [14]. However, decay chains with widely separated half-lives lead to stiff systems of equations. Stiffness makes standard explicit solvers inefficient or unstable. As a result, specialized stiff ODE solvers are often required. Combining analytical insight with numerical methods improves robustness.

1.3 Role of Numerical and Stochastic Simulation

Numerical simulation plays a crucial role when analytical solutions become difficult to evaluate safely or efficiently. Modern ODE solvers can accurately propagate decay systems across long time spans while controlling numerical error [15]. Non-stiff solvers are efficient for short chains with similar half-lives, while stiff solvers handle multi-scale decay behavior more reliably. In addition to deterministic approaches, Monte Carlo simulation models decay as a probabilistic event process. This stochastic method simulates individual decay events using random sampling based on decay probabilities. Monte Carlo models provide physical intuition and serve as an independent validation tool [16]. They are especially useful for demonstrating statistical fluctuations and discrete behavior. Comparing deterministic and stochastic outputs strengthens confidence in results. Hybrid simulation frameworks benefit from both viewpoints. Such integration is increasingly common in scientific computing.

1.4 Need for Stability, Sensitivity, and Visualization

A major requirement in decay chain simulation is numerical stability and physical consistency of computed populations. Negative nuclide counts and overflow errors must be prevented through safe formulations and solver constraints. Sensitivity analysis is also important because uncertainty in decay constants can significantly change long-term predictions [17]. By perturbing parameters and measuring outcome variation, model robustness can be quantified. Energy release and activity calculations extend decay modeling into practical engineering metrics. Visualization of populations, activities, and energy rates improves interpretability and communication of results. Graphical comparisons between analytical, numerical, and Monte Carlo solutions reveal solver behavior clearly [18]. Sensitivity maps further highlight dominant decay parameters. A well-designed simulator therefore includes safety enforcement, multi-method solvers, and rich visualization. These features together support research, teaching, and applied nuclear analysis.

1.5 Computational Implementation Strategy

An effective radioactive decay chain simulator requires a structured computational design that separates physical parameters, numerical solvers, and analysis modules. The decay constants, half-lives, and decay energies should be defined in configurable data structures to allow flexible chain construction [19]. A transition matrix representation enables compact and scalable formulation of the coupled decay equations. Time discretization must be selected carefully to balance resolution and computational cost. Adaptive time solvers further improve efficiency by adjusting step size automatically. Non-negativity constraints should be enforced to maintain physical realism of nuclide populations. Modular solver blocks allow switching between analytical Bateman evaluation and numerical ODE integration. A Monte Carlo engine should be implemented with statistically correct decay probabilities per time step. Reproducibility can be ensured by controlled random number seeding. Such a modular architecture supports reuse and future model extensions.

1.6 Analytical vs Numerical Solver Comparison

Comparing analytical and numerical solutions provides valuable insight into solver performance and model reliability. The Bateman analytical solution offers direct evaluation and serves as a reference baseline for linear chains. However, numerical cancellation and denominator instability can occur when decay constants are similar [20]. Numerical ODE solvers avoid some of these issues by integrating the system dynamically. Explicit solvers such as Runge–Kutta methods perform well for non-stiff cases with moderate time scales. Implicit stiff solvers are better suited for chains with very short and very long half-lives combined. Solver comparison helps identify accuracy, speed, and stability tradeoffs. Cross-validation between methods increases confidence in computed populations. Differences can also reveal discretization sensitivity. Therefore, multi-solver comparison is an important verification step.

1.7 Monte Carlo Validation and Physical Insight

Monte Carlo decay simulation adds a discrete, event-based perspective to decay chain modeling. Instead of solving equations, this method samples decay events using binomial or Poisson statistics. Each nuclide population is updated based on probabilistic decay counts within each time interval [21]. This approach naturally preserves integer particle counts and non-negative behavior. Statistical fluctuations observed in Monte Carlo runs reflect real physical randomness. Ensemble averaging across runs converges toward analytical predictions. Monte Carlo methods are particularly useful for educational demonstration of decay randomness. They also help validate deterministic solver outputs. Differences highlight stochastic variance versus numerical error. Thus, Monte Carlo simulation strengthens physical interpretation and model credibility.

1.8 Applications and Practical Impact

Reliable radioactive decay chain simulation has wide practical importance across science and engineering domains. In nuclear reactor analysis, decay chains determine fuel burnup products and residual heat generation. In radiation protection, activity prediction supports exposure and shielding calculations. Medical physics uses decay modeling for dose planning and isotope production timing. Environmental studies apply decay chains to track radionuclide transport and persistence [22]. Space missions consider radioactive decay for power systems and detector calibration. Geological dating relies directly on long decay chains and daughter isotope ratios. Educational environments use simulation tools to teach nuclear processes interactively. Sensitivity analysis supports uncertainty quantification in safety assessments. Energy release modeling links decay behavior to thermal effects. Therefore, advanced decay chain simulators provide both academic and industrial value.

  1. Problem Statement

Modeling radioactive decay chains accurately remains a challenging problem due to the presence of multiple coupled nuclides with widely different half-lives and decay constants. These systems produce stiff differential equations that often lead to numerical instability and solver inefficiency when handled with basic methods. Analytical Bateman solutions, while exact in theory, can suffer from numerical cancellation and instability in practical computation. Many existing simulation approaches also fail to enforce physical constraints such as non-negative nuclide populations. Deterministic solvers alone may not capture the stochastic nature of radioactive decay events. Furthermore, integrated evaluation of activity, energy release, and parameter sensitivity is often missing in simple models. There is a need for a unified framework that combines analytical, numerical, and stochastic techniques for cross-verification. Such a framework should remain stable across long simulation times and stiff decay chains. It should also provide clear visualization and comparative solver analysis. Therefore, developing a numerically safe and comprehensive radioactive decay chain simulator is the core problem addressed in this work.

You can download the Project files here: Download files now. (You must be logged in).

  1. Mathematical Approach

The mathematical approach for modeling radioactive decay chains is based on coupled first-order differential equations, where the rate of change of each nuclide population depends on its decay and production from its parent. Analytical solutions are obtained using the Bateman equations, which express each nuclide population as a sum of exponential decay terms. For numerical treatment, the system is represented in matrix form and integrated using ODE solvers, allowing both stiff and non-stiff handling. Monte Carlo methods simulate decay probabilistically by sampling decay events at each time step. This combined approach ensures physically consistent, numerically stable, and accurate modeling of complex decay chains. The radioactive decay chain is mathematically described by a set of coupled first-order differential equations:

Analytical solutions are obtained using the Bateman equations:

For numerical modeling, the system is represented in matrix form:

And solved using stiff or non-stiff ODE solvers, while Monte Carlo simulations model decay events probabilistically as binomial processes. The equations describe how each nuclide in a radioactive decay chain changes over time. The first nuclide decreases only through its own decay, while each subsequent nuclide gains population from the decay of its parent and loses population through its own decay. The last nuclide in the chain accumulates all contributions from the previous nuclide without further decay if it is stable. The Bateman equations provide an explicit formula for each nuclide population as a weighted sum of exponential decay terms based on all upstream nuclides. The product and sum terms account for the chain of decays leading to the current nuclide, ensuring the correct contribution from each parent. This formulation guarantees that the populations are physically consistent and non-negative. The matrix form is a compact representation of the entire system, enabling numerical integration with ODE solvers. Stiff solvers handle cases with very different half-lives efficiently, while non-stiff solvers are suitable for simpler chains. Monte Carlo simulations interpret the same equations probabilistically, simulating individual decay events. Together, these approaches provide a complete, accurate, and flexible modeling framework for complex decay chains.

  1. Methodology

The methodology of this study integrates analytical, numerical, and stochastic approaches to simulate radioactive decay chains comprehensively. First, the decay chain is defined by specifying nuclide names, half-lives, decay constants, initial populations, and energy released per decay.

Table 2: Nuclide Parameters

NuclideHalf-life (s)Decay Constant (1/s)Energy per Decay (MeV)
N15.000e+031.386e-045.200
N22.000e+043.466e-051.800
N31.000e+056.931e-060.900
N44.000e+051.733e-060.100
Stableinf0.000e+000.000

Analytical modeling is performed using the Bateman equations, which calculate each nuclide population as a sum of exponential terms weighted by contributions from all parent nuclides. To ensure physical realism, populations are constrained to be non-negative [23]. A matrix representation of the decay chain is constructed, where diagonal elements represent decay of each nuclide and off-diagonal elements capture production from parent nuclides. This formulation allows the system of differential equations to be solved numerically using both non-stiff and stiff ODE solvers, providing flexibility across chains with widely varying half-lives. Solver options are carefully configured with strict relative and absolute tolerances and non-negativity enforcement. Monte Carlo simulation is also implemented to model decay events probabilistically, using binomial sampling for each time step [24]. This stochastic method validates the deterministic solutions and captures natural statistical fluctuations. Activity and total decay rates are computed by multiplying nuclide populations with their respective decay constants. Energy release is calculated by weighting activities with decay energies, providing insight into the thermal and radiological impact of the chain. Sensitivity analysis is performed by perturbing decay constants and observing changes in final nuclide populations, identifying dominant parameters. A fine time grid is chosen to balance computational efficiency and resolution of decay dynamics. Multiple visualizations, including population curves, activity plots, log-scale populations, energy release graphs, solver comparisons, Monte Carlo vs analytical plots, and sensitivity maps, are generated to interpret results effectively [25]. All numerical calculations are safeguarded against negative populations or numerical instability. Random number seeding ensures reproducibility of Monte Carlo simulations. The methodology emphasizes modularity, allowing future extensions to more complex or branched decay chains. This integrated approach provides a robust, accurate, and flexible framework for analyzing radioactive decay chains in research, teaching, and applied nuclear science contexts.

  1. Design Matlab Simulation and Analysis

This simulation models a multi-stage radioactive decay chain using a unified framework that combines analytical, numerical, and stochastic methods to ensure accuracy and physical consistency. The model begins by defining a sequence of nuclides, their half-lives, decay constants, decay energies, and initial populations, forming a linear parent-to-daughter chain ending in a stable isotope. A high-resolution time grid is created to observe both short-term and long-term decay behavior. The analytical part uses a numerically safe Bateman formulation to compute closed-form population solutions while enforcing non-negative values for physical validity. A transition matrix is then constructed to represent decay losses and daughter production in coupled differential equation form. This matrix system is solved using both a non-stiff solver and a stiff solver to compare numerical performance across time scales. Strict solver tolerances and non-negativity constraints are applied to prevent instability and unphysical results. Nuclide activity is computed from decay constants and populations to measure decay rates over time. Energy release rates are also calculated by weighting activities with decay energies, giving insight into total emitted energy. A Monte Carlo module simulates decay events probabilistically at each time step using binomial sampling, capturing stochastic fluctuations. The stochastic results are stored and later compared with analytical predictions for validation. Sensitivity analysis is performed by slightly perturbing each decay constant and recomputing the final populations.

Table 3: Final Populations

NuclideFinal Population
N10.00
N20.00
N35113.49
N4348386.83
Stable646499.68

This reveals which decay parameters most strongly influence system outcomes. Multiple visualizations are generated, including population curves, activity plots, logarithmic population trends, and energy release graphs. Additional plots compare Monte Carlo and analytical results for selected nuclides. Solver comparison graphs illustrate differences between stiff and non-stiff integration. A sensitivity heat map highlights parameter impact across nuclides. All computations include safeguards to prevent negative populations or numerical warnings. The overall simulation provides a robust and extensible platform for studying radioactive decay chain dynamics.

Figure 2: Time evolution of nuclide populations in the radioactive decay chain using the Bateman analytical solution.

You can download the Project files here: Download files now. (You must be logged in).

This figure shows how the population of each nuclide changes over time based on the analytical Bateman solution. The parent nuclide decreases exponentially as it undergoes radioactive decay. Each daughter nuclide first increases due to production from its parent and then decreases as it also decays. Intermediate nuclides exhibit a rise-and-fall behavior typical of decay chains. The stable nuclide continuously accumulates because it has no further decay. The curves demonstrate conservation of total particles across the chain. The smooth behavior confirms numerical stability of the analytical implementation. Non-negative enforcement prevents unphysical values. Population crossover points indicate dominance shifts between nuclides. Overall, the plot validates correct chain dynamics.

Figure 3: Individual nuclide activities versus time showing decay rate variation across the chain.

This figure presents the activity of each nuclide, defined as decay rate over time. Activity depends on both population size and decay constant. Fast-decaying nuclides show sharp early peaks in activity. Longer-lived nuclides produce broader, delayed activity curves. Daughter activities peak after their populations build up. The activity curves reveal when each nuclide contributes most to radiation output. The stable nuclide shows zero activity as expected. The separation of curves illustrates half-life differences clearly. The results are consistent with radioactive decay theory. This plot is useful for radiation intensity assessment.

Figure 4: Log-scale plot of nuclide populations highlighting multi-scale radioactive decay behavior.

This figure displays nuclide populations on a logarithmic scale to reveal multi-scale behavior. Log scaling makes both large and small populations visible simultaneously. Early rapid decay and long-term slow decay appear clearly together. Straight-line segments indicate exponential decay regions. Small daughter populations are no longer hidden near zero. The safe lower bound avoids log-of-zero numerical issues. Differences in decay slopes reflect half-life variation. Transitional growth phases are easier to observe. The plot is especially helpful for stiff decay chains. It confirms stable long-time numerical behavior.

Figure 5: Total energy release rate over time from the complete radioactive decay chain.

This figure shows the total energy release rate from all decays combined. Energy is computed by weighting each activity with its decay energy. Early time energy is dominated by short half-life nuclides. Later energy output shifts to longer-lived daughters. The curve decreases gradually as total activity declines. Small bumps may appear when daughter activity peaks. The plot links decay kinetics with energetic impact. It is useful for thermal and shielding analysis. The smooth curve indicates consistent activity computation. It summarizes the chain’s radiological power output.

Figure 6: Comparison of Monte Carlo stochastic simulation and Bateman analytical results for the second nuclide population.

You can download the Project files here: Download files now. (You must be logged in).

This figure compares stochastic Monte Carlo results with analytical Bateman predictions. The analytical curve is smooth and deterministic. The Monte Carlo curve shows small statistical fluctuations. Both follow the same overall trend over time. Agreement between curves validates the simulation methods. Deviations reflect natural randomness of decay events. Larger populations reduce stochastic noise. The comparison confirms correctness of probability modeling. It demonstrates convergence of stochastic to analytical behavior. This plot supports cross-method verification.

Figure 7: Sensitivity map showing the effect of decay constant variations on final nuclide populations.

This figure presents a sensitivity heat map of final populations to decay constants. Each column represents a perturbed decay parameter. Each row shows the response of a nuclide population. Higher intensity indicates stronger influence. Parent decay constants affect many downstream nuclides. Later-stage constants mainly affect nearby daughters. The map highlights dominant parameters in the chain. It supports uncertainty and robustness analysis. Engineers can identify critical decay data from this view. It guides parameter refinement priorities.

Figure 8: Solver performance comparison between non-stiff and stiff ODE methods for parent nuclide decay.

This figure compares non-stiff and stiff ODE solver results for the parent nuclide. Both solvers produce nearly identical population curves. Close overlap indicates numerical agreement. Minor differences may appear at early rapid decay phases. The stiff solver remains stable for wide time-scale differences. The non-stiff solver performs well in this configuration. The comparison demonstrates solver reliability. It also justifies using stiff solvers for harder chains. Agreement validates matrix formulation correctness. This plot verifies numerical integration quality.

Figure 9: Total radioactive activity evolution for the full decay chain system.

This figure shows the sum of activities from all nuclides over time. It represents total decay rate of the system. The curve decreases as radioactive material is depleted. Early times show higher combined activity. Daughter buildup can slightly delay the decline. The smooth trend confirms consistent activity aggregation. It reflects the overall radiation emission profile. This metric is important for exposure estimation. The result matches expected decay-chain behavior. It provides a single global performance indicator.

  1. Results and Discussion

The simulation results demonstrate that the integrated analytical, numerical, and stochastic framework accurately captures the full dynamics of the radioactive decay chain across multiple time scales.

Table 4: Population Snapshots

Time (s)N1N2N3N4StableTotal Activity
01000000.00.00.00.00.01.386e+02
1597870.05113.6427853.8499805.667227.14.009e+00
3198400.019.4143023.5625704.6231252.62.076e+00
4798930.00.147122.2555811.5397066.21.290e+00
6399470.00.015522.9448071.8536405.38.840e-01
8000000.00.05113.5348386.8646499.76.392e-01

The Bateman analytical solution produced smooth and physically consistent population curves, confirming correct implementation of the closed-form decay relations. Numerical ODE solutions from both non-stiff and stiff solvers closely matched the analytical results, with the stiff solver showing superior robustness for long time spans and widely separated half-lives [26]. All computed populations remained non-negative, validating the enforced physical constraints and solver settings. Intermediate nuclides exhibited the expected growth-and-decay behavior, while the final stable nuclide accumulated monotonically. Activity plots showed clear peak shifting from short-lived parents to longer-lived daughters, illustrating temporal redistribution of radiation output [27]. The total activity curve decreased overall but reflected daughter buildup effects, which is consistent with decay chain theory. Energy release calculations followed activity trends and highlighted the dominant contribution of early high-energy decays. Log-scale population plots confirmed stable numerical behavior even when populations spanned several orders of magnitude. Monte Carlo simulation results aligned well with analytical predictions, with only small statistical fluctuations, supporting probabilistic model correctness. Solver comparison indicated negligible deviation between methods under tight tolerances, strengthening confidence in numerical integration accuracy [28]. Sensitivity analysis revealed that early-chain decay constants have the largest influence on final isotope distributions. Later-stage parameters mainly affected nearby descendants, showing localized sensitivity structure. The sensitivity map provided a clear visual ranking of parameter importance. Cross-validation between analytical, numerical, and stochastic approaches reduced the likelihood of hidden modeling errors. Visualization outputs enabled intuitive interpretation of decay timing and dominance transitions. The framework handled stiff behavior without instability or negative populations. Computational performance remained efficient for the selected resolution. Overall, the results confirm that the proposed simulator is accurate, stable, and suitable for advanced radioactive decay chain analysis.

  1. Conclusion

In conclusion, the presented radioactive decay chain simulator provides a robust and numerically stable framework for modeling multi-stage nuclear decay processes. By combining the Bateman analytical solution with both stiff and non-stiff ODE solvers, the approach ensures high accuracy across different time scales and decay regimes. The inclusion of Monte Carlo simulation adds a stochastic perspective and validates deterministic predictions through probabilistic modeling. Physical constraints such as non-negative populations are consistently enforced, improving realism and reliability [29]. The simulator successfully computes nuclide populations, activities, and energy release profiles in an integrated manner. Sensitivity analysis further identifies the most influential decay parameters, supporting uncertainty evaluation. Comparative solver results demonstrate strong agreement and numerical robustness. Visualization outputs enhance interpretability and educational value [30]. The modular structure allows easy extension to more complex decay networks. Overall, the framework serves as an effective tool for research, teaching, and applied nuclear system analysis.

  1. References

[1] H. Bateman, The solution of a system of differential equations occurring in the theory of radioactive transformations, Proceedings of the Cambridge Philosophical Society, 1910.

[2] G. F. Knoll, Radiation Detection and Measurement, Wiley, 4th edition, 2010.

[3] J. R. Lamarsh and A. J. Baratta, Introduction to Nuclear Engineering, Prentice Hall, 3rd edition, 2001.

[4] D. J. Griffiths, Introduction to Elementary Particles, Wiley, 2nd edition, 2008.

[5] N. Tsoulfanidis and S. Landsberger, Measurement and Detection of Radiation, CRC Press, 4th edition, 2015.

[6] I. Lux and L. Koblinger, Monte Carlo Particle Transport Methods, CRC Press, 1991.

[7] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Springer, 1996.

[8] L. F. Shampine, Numerical Solution of Ordinary Differential Equations, Chapman and Hall, 1994.

[9] C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, 1971.

[10] W. E. Boyce and R. C. DiPrima, Elementary Differential Equations and Boundary Value Problems, Wiley, 10th edition, 2012.

[11] J. D. Jackson, Classical Electrodynamics, Wiley, 3rd edition, 1998.

[12] S. Wolfram, The Mathematica Book, Wolfram Media, 5th edition, 2003.

[13] MATLAB Documentation, Ordinary Differential Equation Solvers, MathWorks.

[14] MATLAB Documentation, Stiff ODE Solvers and Numerical Stability, MathWorks.

[15] R. L. Burden and J. D. Faires, Numerical Analysis, Brooks Cole, 9th edition, 2011.

[16] K. E. Atkinson, An Introduction to Numerical Analysis, Wiley, 2nd edition, 1989.

[17] J. M. Thijssen, Computational Physics, Cambridge University Press, 2nd edition, 2007.

[18] W. H. Press et al., Numerical Recipes in Scientific Computing, Cambridge University Press, 3rd edition, 2007.

[19] P. Arbenz, Introduction to Parallel Computing, Oxford University Press, 2015.

[20] J. Spanier and E. M. Gelbard, Monte Carlo Principles and Neutron Transport Problems, Dover, 2008.

[21] D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, 2nd edition, 2004.

[22] M. E. J. Newman, Computational Physics, CreateSpace, 2013.

[23] A. Quarteroni, R. Sacco, and F. Saleri, Numerical Mathematics, Springer, 2nd edition, 2007.

[24] J. W. Demmel, Applied Numerical Linear Algebra, SIAM, 1997.

[25] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, 2nd edition, 2003.

[26] R. G. Edwards, Numerical Methods for Physics Applications, Wiley, 2018.

[27] P. J. Mohr, D. B. Newell, and B. N. Taylor, CODATA recommended values of the fundamental physical constants, Reviews of Modern Physics, 2016.

[28] International Atomic Energy Agency, Nuclear Data Services Handbook, IAEA.

[29] National Nuclear Data Center, Chart of Nuclides and Decay Data Resources.

[30] J. J. Duderstadt and L. J. Hamilton, Nuclear Reactor Analysis, Wiley, 1976.

You can download the Project files here: Download files now. (You must be logged in).

Related Articles

Responses

Your email address will not be published. Required fields are marked *

L ading...