{ "cells": [ { "cell_type": "markdown", "id": "18fd7414", "metadata": {}, "source": [ "# Lesson 5 - Monte Carlo First" ] }, { "cell_type": "markdown", "id": "a749ae9f", "metadata": {}, "source": [ "## Objectives\n", "\n", "After this lesson, you should be able to do the following:\n", "\n", "- Explain the difference between a **model**, a **simulation**, and an **observable**.\n", "- Distinguish a **deterministic** view from a **stochastic** view of the same physical system.\n", "- Use a random number as a computational coin for a yes/no event.\n", "- Describe the basic Monte Carlo pattern: one history, one tally, many histories, one estimate.\n", "- Read and write the logic for a simple neutron source-detector Monte Carlo simulation.\n", "- Separate parameters that change the physical model from parameters that change the sampling noise." ] }, { "cell_type": "markdown", "id": "400c4581", "metadata": {}, "source": [ "## Why Monte Carlo first?\n", "\n", "When I took on this course, I knew we'd need to introduce probability and statistics to support later courses. Rather than treat these as a separate language (like many math courses, in my experience), we'll consider an application for which probability and statistics offer an enabling toolbox.\n", "\n", "But that starts with the history of one neutron. \n", "\n", "This neutron is born at a source. It moves through a material. It might bounce off a nucleus. It might be absorbed. It might escape the volume of interest. It might reach a detector. One neutron has one life. A detector response is not the life of one neutron. It is a statement about many possible neutron lives.\n", "\n", "That is the core reason Monte Carlo simulation is fundamental to nuclear engineering. Particle transport is naturally a story of individual histories. We can often write equations for the average behavior, and those equations are important. But the event-by-event view is also powerful because it mirrors the physical narrative:\n", "\n", "```text\n", "source particle\n", " -> random flight or step\n", " -> possible interaction\n", " -> possible termination\n", " -> tally\n", " -> repeat many times\n", " -> estimate response\n", "```\n", "\n", "In later courses on shielding and nuclear reactor physics, those histories will become very, very complex. A neutron will have an energy, a direction, a position, material-dependent cross sections, scattering kernels, and detector response functions. For now, we will build the smallest version of the idea.\n", "\n", "This lesson is called **Monte Carlo first** because we will use randomness before formalizing all of probability. We will use the language of events, tallies, histories, and estimates. In the next lessons, we will give more formal names to the probability objects that are already doing work here." ] }, { "cell_type": "markdown", "id": "5626aa4f", "metadata": {}, "source": [ "## Models, simulations, and observables\n", "\n", "A **model** is a representation of a physical system. It is the collection of assumptions, balance equations, probabilistic rules, parameters, data, and numerical approximations used to describe what matters and how those pieces relate.\n", "\n", "A model is not the real system. It is a chosen description of the real system.\n", "\n", "A **simulation** is the execution of a model. It is what happens when the model is run to produce predictions, histories, fields, tallies, or other outputs.\n", "\n", "The model is the representation. The simulation is the run.\n", "\n", "Before building either one, we should identify the **observable**. The observable is the quantity we want to predict, estimate, or understand. Examples include transmitted intensity, remaining nuclei, detector counts, neutron flux, absorbed dose, leakage probability, and probability of detection.\n", "\n", "This distinction matters because the same physical system can be modeled in more than one way. A model is a choice about what is kept, what is averaged, what is ignored, and what is sampled.\n", "\n", "For example, suppose a beam travels through an attenuating slab. If the observable is only the transmitted fraction, then one model may describe a smooth intensity field. Another model may follow individual particles. Both can aim at the same observable.\n", "\n", "The observable keeps us honest. Without it, we can add details forever and still not know what question we are answering." ] }, { "cell_type": "markdown", "id": "d2e6a85d", "metadata": {}, "source": [ "## Deterministic and stochastic views\n", "\n", "A **deterministic** model has no random branching. Once the inputs are fixed, the output is fixed. If you run the same deterministic model twice with the same parameters, initial conditions, boundary conditions, and numerical settings, you get the same result.\n", "\n", "A deterministic model often describes an averaged quantity. Things like temperature, pressure, and concentration are commonly treated this way. Many microscopic events have already been averaged into the input parameters.\n", "\n", "A **stochastic** model includes random events. The same inputs do not define one path through the model. They define a collection of possible paths and how likely each path is.\n", "\n", "A Monte Carlo simulation is a stochastic simulation that estimates an observable by sampling many random histories. One history is one possible realization of the model. The answer comes from many histories summarized together.\n", "\n", "Despite their differences, deterministic and Monte Carlo simulation are often two routes to the same physical quantity.\n", "\n", "The deterministic route averages first, then solves for a smooth response. The Monte Carlo route samples individual events first, then averages afterward." ] }, { "cell_type": "markdown", "id": "9792d842", "metadata": {}, "source": [ "## A first physical example: attenuation through a slab" ] }, { "cell_type": "markdown", "id": "de6c3006", "metadata": {}, "source": [ "### The physical system\n", "\n", "A narrow beam enters an attenuating material of thickness $L$. The material might be tinted glass, fog, murky water, shielding, or any medium that removes some of the beam as it travels.\n", "\n", "The observable is the **transmitted fraction**:\n", "\n", "$$\n", "T = \\frac{\\text{amount leaving the far side}}{\\text{amount entering the slab}}.\n", "$$\n", "\n", "We will describe the same system in two ways.\n", "\n", "First, we will use a deterministic continuum model. Then, we will use a stochastic particle model." ] }, { "cell_type": "markdown", "id": "71a2d3d6", "metadata": {}, "source": [ "### Deterministic continuum model\n", "\n", "Let $I(x)$ be the beam intensity after traveling distance $x$ into the material. Suppose the material removes a fixed fraction of the remaining intensity per unit length. Then the intensity satisfies\n", "\n", "$$\n", "\\frac{dI}{dx} = -\\mu I,\n", "$$\n", "\n", "where $\\mu$ is the attenuation coefficient; such coefficients are almost always evaluated based on a combination of measurements and modeling. \n", "\n", "The minus sign says the intensity decreases as the beam travels. The factor $I$ says the amount removed over a small distance is proportional to how much remains.\n", "\n", "With $I(0)=I_0$, the solution is\n", "\n", "$$\n", "I(x)=I_0 e^{-\\mu x}.\n", "$$\n", "\n", "At the far side of the slab,\n", "\n", "$$\n", "I(L)=I_0 e^{-\\mu L},\n", "$$\n", "\n", "so the transmitted fraction is\n", "\n", "$$\n", "T = \\frac{I(L)}{I_0}=e^{-\\mu L}.\n", "$$\n", "\n", "This model does not follow individual particles. It treats intensity as a smooth, continuous quantity. It is a continuum model and gives the average response directly." ] }, { "cell_type": "markdown", "id": "889e4c26", "metadata": {}, "source": [ "### Stochastic particle model\n", "\n", "Now imagine launching identical particles into the slab one at a time. Each particle has a history. It either reaches the far side or it is removed inside the material.\n", "\n", "For one particle, define a yes/no outcome:\n", "\n", "$$\n", "Y =\n", "\\begin{cases}\n", "1, & \\text{if the particle reaches } x=L,\\\\\n", "0, & \\text{if the particle is removed before reaching } x=L.\n", "\\end{cases}\n", "$$\n", "\n", "A single particle does not give the transmitted fraction. A single particle gives only one yes/no result.\n", "\n", "To connect this event-by-event picture to the continuum model, divide the slab into many thin layers of width $\\Delta x$. Suppose a particle entering one layer has probability approximately $\\mu\\Delta x$ of being removed in that layer and probability $1-\\mu\\Delta x$ of surviving to the next layer.\n", "\n", "If the slab contains\n", "\n", "$$\n", "M = \\frac{L}{\\Delta x}\n", "$$\n", "\n", "layers, then the probability of surviving all layers is approximately\n", "\n", "$$\n", "(1-\\mu\\Delta x)^M.\n", "$$\n", "\n", "As the layers become thinner,\n", "\n", "$$\n", "(1-\\mu\\Delta x)^{L/\\Delta x} \\rightarrow e^{-\\mu L}.\n", "$$\n", "\n", "Thus the stochastic particle model and the deterministic attenuation equation predict the same transmitted fraction in the thin-layer limit.\n", "\n", "The deterministic model gives $T=e^{-\\mu L}$ directly. The stochastic model estimates that same value by simulating many particle histories.\n", "\n", "For $N$ simulated particles,\n", "\n", "$$\n", "\\hat{T} = \\frac{1}{N}\\sum_{i=1}^{N}Y_i\n", "= \\frac{\\text{number of transmitted particles}}{\\text{number of launched particles}}.\n", "$$\n", "\n", "The symbol $T$ denotes the model quantity. The symbol $\\hat{T}$ denotes the estimate from one finite simulation.\n", "\n", "That distinction is central to Monte Carlo methods:\n", "\n", " - $T$ is what the model predicts in the ideal infinite-sampling limit. \n", " - $\\hat{T}$ is what we observe with a finite number of histories." ] }, { "cell_type": "markdown", "id": "fccde8bd", "metadata": {}, "source": [ "## Random numbers as computational coins\n", "\n", "A computer can generate pseudorandom numbers between 0 and 1. In Python, we often get them from a NumPy random number generator:" ] }, { "cell_type": "code", "execution_count": 1, "id": "21394053", "metadata": { "execution": { "iopub.execute_input": "2026-06-08T15:52:13.651887Z", "iopub.status.busy": "2026-06-08T15:52:13.651613Z", "iopub.status.idle": "2026-06-08T15:52:13.789960Z", "shell.execute_reply": "2026-06-08T15:52:13.789552Z" } }, "outputs": [], "source": [ "import numpy as np\n", "\n", "rng = np.random.default_rng()\n", "U = rng.random()" ] }, { "cell_type": "markdown", "id": "eb5c5c25", "metadata": {}, "source": [ "The value $U$ behaves like a random draw from the interval $[0,1)$. We might think of $U$ as a computational, $n$-sided die whose possible face values $v = i/n, i = 0, 1, 2, \\ldots, n-1$ uniformly fill that interval as $n$ grows very, very large.\n", "\n", "Suppose an event should occur with probability $p$, where $0\\le p\\le 1$. Draw $U$ between 0 and 1 and use the rule\n", "\n", "$$\n", "U < p.\n", "$$\n", "\n", "If $U= n_steps:\n", " return 1\n", " if x < 0:\n", " return 0\n", " # If still inside, sample absorption.\n", " if rng.random() < p_absorb:\n", " return 0\n", " # Otherwise the neutron scattered and continues.\n", "\n", " # It is good practice to report unexpected issues!\n", " print(\"Warning: history did not terminate before max_steps!\")" ] }, { "cell_type": "markdown", "id": "5fdf51ef", "metadata": {}, "source": [ "The function returns a single score. It does not return the whole path. That is because the observable is detector response. If the observable were path length, number of collisions, or energy deposited, the tally would be different.\n", "\n", "Now repeat many histories:" ] }, { "cell_type": "code", "execution_count": null, "id": "cf41da70", "metadata": { "execution": { "iopub.execute_input": "2026-06-08T15:52:13.802698Z", "iopub.status.busy": "2026-06-08T15:52:13.802512Z", "iopub.status.idle": "2026-06-08T15:52:13.805236Z", "shell.execute_reply": "2026-06-08T15:52:13.804806Z" } }, "outputs": [], "source": [ "def run_simulation(n_histories, n_steps, p_absorb, seed=None):\n", " \"\"\"Estimate the detector response from many neutron histories.\"\"\"\n", " rng = np.random.default_rng(seed)\n", "\n", " detected = 0\n", "\n", " for _ in range(n_histories):\n", " detected += run_one_history(\n", " n_steps=n_steps,\n", " p_absorb=p_absorb,\n", " rng=rng,\n", " )\n", "\n", " return detected / n_histories" ] }, { "cell_type": "markdown", "id": "fbffdd0a", "metadata": {}, "source": [ "An example run might look like this:" ] }, { "cell_type": "code", "execution_count": 4, "id": "48d6a44d", "metadata": { "execution": { "iopub.execute_input": "2026-06-08T15:52:13.806877Z", "iopub.status.busy": "2026-06-08T15:52:13.806505Z", "iopub.status.idle": "2026-06-08T15:52:13.830826Z", "shell.execute_reply": "2026-06-08T15:52:13.830368Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0059\n" ] } ], "source": [ "response = run_simulation(\n", " n_histories=10_000,\n", " n_steps=10,\n", " p_absorb=0.1,\n", " seed=12345,\n", ")\n", "\n", "print(response)" ] }, { "cell_type": "markdown", "id": "f5bc00e4", "metadata": {}, "source": [ "If you change the seed, the answer may change slightly. That is not a bug. It is sampling noise.\n", "\n", "If you increase `n_histories`, the estimate will usually wiggle less from run to run. That does not change the physical model. It changes the quality of the estimate." ] }, { "cell_type": "markdown", "id": "5fe08de7", "metadata": {}, "source": [ "## Reading the output\n", "\n", "The output of the simulation is an estimate:\n", "\n", "$$\n", "\\hat{p}_d = \\frac{\\text{number detected}}{\\text{number simulated}}.\n", "$$\n", "\n", "It is tempting to say, \"the detector probability is $\\hat{p}_d$.\" A sharper statement is: \"For this model and this finite simulation run, the estimated detector response was $\\hat{p}_d$.\"\n", "\n", "The actual model probability $p_d$ is the ideal quantity we are trying to estimate. We do not know it exactly just because we ran a finite simulation.\n", "\n", "This is not a philosophical detail. It affects how we compare simulations.\n", "\n", "Suppose two simulation runs give 0.041 and 0.044. Is the second model really more detectable? Or did two finite Monte Carlo estimates simply wiggle? We cannot answer that without thinking about sampling uncertainty. That topic belongs to the later statistics lesson, but the habit starts here." ] }, { "cell_type": "markdown", "id": "6761a5a8", "metadata": {}, "source": [ "## What changes the model? What changes the run?\n", "\n", "Monte Carlo results can change for different reasons. We should separate them." ] }, { "cell_type": "markdown", "id": "93eeac7e", "metadata": {}, "source": [ "### Model parameters\n", "\n", "Model parameters define the physical or mathematical system being represented.\n", "\n", "In the attenuation example, $\\mu$ and $L$ are model parameters.\n", "\n", "In the neutron toy model, model choices include\n", "\n", "- geometry;\n", "- detector boundary location;\n", "- leakage boundary location;\n", "- step size;\n", "- probability of absorption;\n", "- probability of moving left or right;\n", "- boundary rules;\n", "- interaction rules.\n", "\n", "Changing any of these changes the model. It asks a different physical or mathematical question." ] }, { "cell_type": "markdown", "id": "6089e618", "metadata": {}, "source": [ "### Sampling parameters\n", "\n", "Sampling parameters change how the model is executed. The number of histories $N$ and the random seed are sampling parameters.\n", "\n", "Increasing $N$ usually reduces sampling noise. It does not change the underlying model.\n", "\n", "Changing the seed changes the particular random realization. It should not systematically change the model prediction. If changing the seed creates a large persistent shift, either the simulation is too noisy or the code has a problem.\n", "\n", " " ] }, { "cell_type": "markdown", "id": "bc5b53db", "metadata": {}, "source": [ "## What we have already used\n", "\n", "We have already used several probability and statistics ideas, even though we have not fully formalized them.\n", "\n", "- **Events**: detection, absorption, leakage, scattering, and transmission are possible outcomes of a random process.\n", "- **Indicators**: a yes/no tally equals 1 when the event of interest happens and 0 otherwise.\n", "- **Probability rules**: every random decision is controlled by a rule that assigns likelihoods to possible outcomes; later we will call these distributions.\n", "- **Estimates**: a Monte Carlo result approximates a model quantity, such as a probability, average response, detector count rate, flux, or other observable.\n", "- **Sampling noise**: finite numbers of random histories make estimates wiggle from run to run; later we will quantify that wiggle." ] }, { "cell_type": "markdown", "id": "ccd25c94", "metadata": {}, "source": [ "## Deterministic and Monte Carlo views of the same problem\n", "\n", "Return to attenuation through a slab.\n", "\n", "The deterministic model solves\n", "\n", "$$\n", "\\frac{dI}{dx}=-\\mu I\n", "$$\n", "\n", "and obtains\n", "\n", "$$\n", "\\frac{I(L)}{I_0}=e^{-\\mu L}.\n", "$$\n", "\n", "The stochastic model launches many particles and estimates the same transmitted fraction from the number that reach the far side.\n", "\n", "These are not unrelated tricks. They are two ways of arranging the averaging.\n", "\n", "The deterministic view says: average the beam into a smooth intensity and solve for the average response.\n", "\n", "The Monte Carlo view says: simulate individual histories, tally the outcome of each history, then average the tallies.\n", "\n", "The neutron source-detector model follows the same pattern. A deterministic transport method would try to describe the average particle population or flux directly. A Monte Carlo method follows many particles and estimates the detector response from their histories.\n", "\n", "This is why we start with Monte Carlo. It makes the physical story visible." ] }, { "cell_type": "markdown", "id": "fb41f389", "metadata": {}, "source": [ "## Looking ahead\n", "\n", "The next lessons will formalize the objects we have already used.\n", "\n", "In the next lesson, we will name random variables and distributions. We will describe what `rng.random()` is sampling from, and we will define the probability rules behind random decisions.\n", "\n", "After that, we will learn how to transform uniform random numbers into samples from other distributions. That will let us replace fixed neutron steps with sampled flight distances.\n", "\n", "Finally, we will quantify Monte Carlo uncertainty. We will ask how much a finite estimate should wiggle and how to report that wiggle responsibly." ] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.3" }, "livereveal": { "scroll": true } }, "nbformat": 4, "nbformat_minor": 5 }