Radionuclides Calculation

I recently asked Opus 4.6, Gemini 3 Pro, and GPT 5.2 Thinking to compute some radionuclide data. I have $20/month plans with GPT and Gemini and a $100/month plan with Opus, so it's not a fair comparison. But Opus did perform best.

From the Gemini UI, it wasn't clear whether I should choose Thinking or Pro. From a web search, it seems Thinking is the small model with CoT, whereas Pro is the big model. So I chose Pro.

Claude's split-view UI for artifacts is a big win. I'm amazed the other two haven't copied it. Though according to Claude, it doesn't know which artifact you're viewing, so you have to mention it all the time. How hard could it be to make Claude aware of this?

None of the models could produce a complete list of radionuclides with half-lives 4-40 years. So I made the list myself with NuDat 3 and gave it to them. Fun fact: there are only 24.

I asked the models to compute the decay power density of each nuclide at t=0 and put the results in a table. Then I asked them to add a "P10" column showing the average power density over ten years. This is a relatively hard task which requires looking up decay data for hundreds of nuclides and integrating the decay energy over trees of nuclides as they change over time. I suggested the models write a Monte Carlo simulation.

Gemini performed worst. It got incorrect results for the P10 calculation (I didn't confront it about this). And as usual: would only export tables to Google Sheets, specified units in every cell in the table (rather than in column headers), and often printed invalid LaTeX in its output.

GPT got generally correct results and shared its Monte Carlo code.

Opus gave the best results and happily produced output in markdown, csv, and xlsx.

As a point of reference, I once gave almost exactly this same task to a nuclear engineer (with a PhD) and he came back a week later with only a partial solution.

When I asked Opus to share its Monte Carlo code, it said:

I should be honest: I never actually wrote a Monte Carlo simulator. You suggested that approach, but I went with deterministic ODE integration of the Bateman equations instead, which is exact (up to numerical tolerance) rather than statistical. For systems with well-characterized branching ratios and no stochastic fluctuations to worry about, the ODE approach is both faster and more accurate than Monte Carlo would be.

The Bateman equations describe how a chain of radioactive species evolves over time. The basic idea is simple: each nuclide is created by its parent's decay and destroyed by its own decay. For a chain A → B → C → ... the number of atoms of species i changes as:

$$\frac{dN_i}{dt} = \lambda_{i-1} N_{i-1} - \lambda_i N_i$$

where $\lambda = \ln 2 / t_{1/2}$ is the decay constant. The first term is the feed-in from the parent, the second is loss from its own decay. For the top of the chain there's no feed-in, so $dN_1/dt = -\lambda_1 N_1$, which just gives pure exponential decay.

Bateman published the analytical solution in 1910 — for a simple linear chain starting from pure parent, each $N_i(t)$ is a sum of exponentials:

$$N_i(t) = N_1(0) \cdot \prod_{j=1}^{i} \lambda_j \cdot \sum_{j=1}^{i} \frac{e^{-\lambda_j t}}{\prod_{k \neq j} (\lambda_k - \lambda_j)}$$

That looks ugly but it's just a superposition of decaying exponentials with coefficients determined by all the decay constants in the chain. It works beautifully for a clean linear chain like Ra-228 → Ac-228 → Th-228 → Ra-224 → ...

Why I didn't use the analytical solution:

Three complications make the closed-form Bateman solution impractical for our problem:

First, branching decays. Bi-212 goes 64% to Po-212 and 36% to Tl-208. Ac-227 goes 99% to Th-227 and 1% to Fr-223. The chain isn't linear — it's a tree. Bateman's formula assumes no branching.

Second, extreme stiffness. In the Ra-228 chain, the decay constants span 15 orders of magnitude — Ra-228 at $3.8 \times 10^{-9}$ s⁻¹ versus Po-212 at $2.3 \times 10^{6}$ s⁻¹. When you try to evaluate the analytical formula, you get catastrophic numerical cancellation — enormous exponentials being subtracted from each other to give a tiny result. This is a classic problem with Bateman's formula and it's why practical nuclear engineering codes almost never use it directly.

Third, I need total energy deposited, not just populations. Each species deposits a different amount of energy per decay, and I need the time integral of power over 10 years. Analytically that means integrating $\lambda_i N_i(t) E_i$ for every species — doable but messy with branches.

What I did instead — the ODE approach:

I wrote the general matrix form:

$$\frac{d\vec{N}}{dt} = (-\text{diag}(\vec{\lambda}) + F) \cdot \text{diag}(\vec{\lambda}) \cdot \vec{N}$$

where $F$ is the "feeding matrix" — $F_{ij}$ is the branching fraction from species $j$ to species $i$. This handles any topology: linear chains, branches, even converging chains (like Es-252 feeding both Bk-248 and Cf-252 which then rejoin at Cm-244).

I added one extra variable to the ODE vector: cumulative energy deposited in MeV. Its derivative is just the instantaneous power:

$$\frac{dE}{dt} = \sum_i \lambda_i N_i \cdot \epsilon_i$$

where $\epsilon_i$ is the energy deposited per decay of species $i$. Then $P_{10} = E(10\text{ yr}) / (10\text{ yr})$ falls right out.

The collapsing trick:

Species with $t_{1/2} < 1$ hour (like Rn-220 at 55.6 seconds, or Po-212 at 0.3 μs) reach secular equilibrium essentially instantaneously relative to the timescales we care about. Tracking them in the ODE would force the solver to take microsecond timesteps for 10 years — computationally insane. So I "collapsed" them: if a short-lived daughter decays to another short-lived daughter, I follow the chain until I reach something with $t_{1/2} \geq 1$ hour, accumulating all the deposited energy along the way. For example, Ra-224's effective decay becomes:

Ra-224 → (Rn-220 → Po-216 instantly) → Pb-212, depositing 5.79 + 6.41 + 6.91 = 19.11 MeV

That's why Ra-224 shows 19.11 MeV per decay in the simulation — it's carrying the collapsed energy of three fast α decays.

The solver:

Scipy's Radau (implicit Runge-Kutta) is designed for stiff systems. Even after collapsing, the remaining decay constants still span several orders of magnitude (hours to decades), so an explicit solver like RK45 would struggle. Radau handles this gracefully with adaptive timesteps — large steps when things are slowly evolving, small steps during rapid transients like daughter ingrowth.

So the whole method is: build the chain, collapse the fast stuff, set up the matrix ODE, throw it at Radau, read off the energy at 10 years. Works for any nuclide in about 0.1 seconds.

Pretty... pretty good!

-Carl

clumma/nuclides_test.md

Select an option

No results found

Select an option

No results found

Radionuclides Calculation

See also: