Calculation of K2: D/R² ~ 1e-9 / (1e-6)² = 1000 per second. Fix: Original calculation was incorrect; recalculation of K2 from 1e-9 m²/s is 1000. This implies a millisecond timescale for diffusion out.
K3 (oxidation reaction): Referencing usual rates (10^8 L/mol·s with activation energy ~15 kcal/mol). K3 from the model must be per second units, so it's the lookup rate constant times the RH concentration.
Approximating RH concentration: From density of organic matter and H atoms per CH2 group, unit conversions from g/cm³ to mol/m³. Must use unit conversions.
Estimated value calculation attempt for K3: With body temperature (310K) and estimated values. First calculation appeared too small (5e-25).
Correction: Saw unit/mole problem in RH concentration, corrected calculation. Recalculation of K3 is 30 per second.
K4 (radical termination): Normal rates are 10^5 L/mol·s for recombination of peroxy radicals. K4 in the simulation must factor in vesicle volume. K4 is standard K divided by reactor volume, possibly incorporating Avogadro's number for units of single molecules.
Determination of K4: 10^5 L/mol·s divided by vesicle volume (4/3 * pi * R³) and Avogadro's number, making unit conversions (L to m³).
Attempt at calculation for K4: Student calculation correction. Recalculated K4 is 4e-5 per second.
Determining maximum simulation time (T_max): Guessing 1000 seconds.
Printing the calculated rates: K = [1e-5, 1000, 30, 4e-5].
First concern with the rates: K2 (1000) is too fast, implying that radicals escape rapidly. The originally guessed K1 (1e-5) is extremely slow. This is an issue.
Reassessing K1: K1 must agree with the rate of diffusion in, which must be comparable with the rate of diffusion out (~1000 per second), multiplied by the average concentration of radicals outside. If the average concentration inside/outside is comparable, K1 must be comparable with K2.
K1 guessing: If average internal concentration is 1 radical, K1 must be ~1000. If average internal concentration is 10 radicals, K1 must be ~10000. Let's try K1 = 1000. K = [1000, 1000, 30, 4e-5].
Initial attempt to execute Matlab code.
Fixing code errors: (parentheses, saving file).
Plotting results: (peroxides count vs. time). First plot appears weird (two points only).
Looking at the trajectory data directly: Two points only were measured. Peroxide count jumps to 307, which is crazy for a single step increment. This points to a bug.
Finding the bug: The step counter and path storage were outside the primary time loop, and they only executed once. Placed step update and storage within the loop.
Second try at executing the code: It's executing now ("doing something").
Fearing nested loops or endless loops: The step limit must be enforced in order not to flood memory.
Advice to use: `break` when max steps is reached. By using a `while` loop condition of `(time < Tmax) & (step < max_steps)`.
Resolving a possible bug: Using the `&&` operator in the while loop condition.
Third try running the code: Too quick, says it completed instantly.
Fourth try: The number of radicals plot vs. time looks better, with jumps up and down as desired. This path exhibits sampling from the stochastic process. The peroxide plot is still problematic.
Simulation Result Analysis
Must run numerous trajectories (e.g., 10,000) since each simulation is a sample from the probability distribution.
Methods of Analyzing Several Trajectories:
- Averaging trajectories: Adding together trajectories and dividing by the number of runs to obtain an average trajectory. Problem: Trajectories have different points in time, requiring a method of compounding them (e.g., interpolation or certain plot functions).
- Histograms: Snapshots at one time point across all trajectories and a plot of the distribution of the number of radicals or peroxides.
- Plotting figures: Such as the amount of time to achieve a given number of peroxides or standard deviation.
- Challenge: Computing desired statistics such as mean and standard deviation from the sampled output.
- Alternative approach: Solving the master equation as a system of ordinary differential equations (ODEs). This gives the probability distribution P(n_rad, n_ro, time) explicitly.
- Data analysis issue: Even for ODE solutions, handling the high-dimensional probability distribution at various time points is crucial. People tend to plot the average number of species versus time (e.g., average n_ro). One should also look at the dispersion (e.g., n_rad squared average, standard deviation).
- Reference: A homework exercise on kinetic Monte Carlo for catalysis, with an extension offer for homework due dates.
- Conclusion: Summary of the lecture.
.png)
.png)