Progressing

docs/src/Teaching/Transport Phenomena/02-Applications-to-mass-and-heat-transfer.md

@@ -68,14 +68,19 @@ The models discussed so far are in fact the macroscopic response to discrete phe
 
 In the most simple of cases, a *non-interacting particle* in a dilute mixture performs a *random walk*; when many such particles are present in one part of the domain, as we will show later, the resulting gradient of concentration that will be established over a sufficiently large amount of time leads to Fick's law. The terminology was created by [Pearson (1905)](https://www.nature.com/articles/072294b0) for which [Rayleigh (1905)](https://www.nature.com/articles/072318a0) had already found a solution many years beforehand; is also appears in the solution of Brownian motion by [Einstein (1905)](https://myweb.rz.uni-augsburg.de/~eckern/adp/history/einstein-papers/1905_17_549-560.pdf) (translated [here](https://www.damtp.cam.ac.uk/user/gold/pdfs/teaching/old_literature/Einstein1905.pdf)) and turbulent diffusion by [Taylor (1922)](https://courses.washington.edu/mengr537/Lecture_Notes/DiffusionContinuousMovements_TaylorDispersion_ProcMathSocLon1921.pdf).
 
-In the derivation of diffusion equation, following the already mentioned *non-interacting particles* hypothesis, one also needs steps to be *independent, identically distributed* (IDD).
+In the derivation of diffusion equation, following the already mentioned *non-interacting particles* hypothesis, one also needs steps to be *independent, identically distributed* (IDD) and limited by a *finite variance* $\sigma^{2}$, *i.e.* $p(x\pm{}K\sigma)\to{}0$ quickly for some finite $K$.
 
 $$
-P_{N+1}(x) = \int{}P_{N}(x-y)p(y)dy
+P_{N+1}(x) = \int{}P_{N}(x-y)p(y)dy = (P_{N} \ast p)(x)
 $$
 
 > *The probability $P_{N+1}(x)$ that a particle is found at $x$ at step $N+1$ is equal to the sum of the probabilities $p(y)$ of particles at a distance $y$ from $x$ execute a movement of such amplitude to reach $x$ times the probability of a particle being at the position $P_{N}(x-y)$.*
 
+$$
+P_{N} = (P_{0} \ast p \ast p \ast \dots p)(x) = (P_{0} \ast p^{\ast{}N})(x)
+$$
+
+
 On the other limit one finds the concentrated mixtures for which particle-particle interactions intervene on the probability distribution of motion and results obtained through *non-equilibrium thermodynamics* lead to a more complex law.
 
 

script/random-walk.jl

@@ -4,16 +4,19 @@ using ProgressMeter
 using Random
 
 # Probility of moving a distance y.
-p(z) = exp(-z / 3)
+p(z, K) = exp(-z / K)
 
 # Create a random number generator with a given seed.
 rng = MersenneTwister(42)
 
 # Size of square domain (matrix) for random walks.
-n = 20
+n = 100
 
 # Number of random walk steps to perform.
-N = 1_000_000
+N = 1_000
+
+# Maximum step size.
+K = 2
 
 # Maximum number of tries for performing a step.
 max_tries = 100
@@ -25,52 +28,67 @@ max_fails = 0.01
 A = zeros(Int64, (n, n))
 
 # Get initial position of *particle*.
-xn, yn = rand(rng, 1:n, 2)
+xn, yn = rand(rng, 0:n-1, 2)
+
+# Store initial position for display later.
 x0, y0 = xn, yn
 
 # Feed particle to its initial position.
-A[xn, yn] += 1
+A[xn+1, yn+1] += 1
 
+# Tracker for number of misses.
 frozen = 0
 
+# Initialize figure.
+# fig, ax, hm = heatmap(A / N; colormap = :thermal)#, colorrange = (0.0, N/n))
+# scatter!(ax, x0, y0; color = :black)
+
+saveas = joinpath(@__DIR__, "random-walk.mp4")
+
 @showprogress for tn ∈ 1:N
+# record(fig, saveas, 1:N) do tn
+    # global xn, yn, A, frozen
+
     n_try = 0
 
     while n_try < max_tries
         n_try += 1
 
-        xm, ym = rand(rng, 1:n, 2)
-
-        # TODO: consider the system wraps around!
-        z = hypot((xn-xm), (yn-ym))
+        # TODO: system wraps around! Not fully working if we still use clamp!
+        Δx, Δy = rand(rng, -K:K, 2)
+        xm = mod(xn + Δx, n)
+        ym = mod(yn + Δy, n)
 
-        bullet = rand()
+        z = hypot(xn-xm, yn-ym)
 
-        threshold = p(z)
+        bullet = rand()
+        threshold = p(z, K)
 
         if bullet < threshold
             xn, yn = xm, ym
-            A[xn, yn] += 1
+            A[xn+1, yn+1] += 1
             break
         end
     end
 
+    # hm[3].val = A / N
+
     if n_try == max_tries
         frozen += 1
     end
 end
 
+# hm[3].val = A / sum(A)
+
 if frozen > max_fails * N
     @error("Missed too many steps: $(frozen) of $(N)")
 end
 
-B = A / sum(A)
-
 fig = with_theme() do
     f = Figure()
     ax = Axis(f[1, 1])
 
-    heatmap!(ax, B; colormap = :thermal)
+    heatmap!(ax, A / N; colormap = :thermal)
     scatter!(ax, x0, y0; color = :black)
     f
 end