From 0042d4e69cf1575dd59d42a2e2579ae597919d41 Mon Sep 17 00:00:00 2001
From: Orestis <orestis.malaspinas@pm.me>
Date: Wed, 13 Mar 2024 23:41:28 +0100
Subject: [PATCH] updates

---
 .../content/1-dimensionless-lbm.md            |  2 +-
 sage/hermite.sage                             | 28 +++++++++++++------
 2 files changed, 21 insertions(+), 9 deletions(-)

diff --git a/lattice-boltzmann.com/content/1-dimensionless-lbm.md b/lattice-boltzmann.com/content/1-dimensionless-lbm.md
index 394f79b..b2cb286 100644
--- a/lattice-boltzmann.com/content/1-dimensionless-lbm.md
+++ b/lattice-boltzmann.com/content/1-dimensionless-lbm.md
@@ -113,7 +113,7 @@ $$
 \partial_t^\ast f^\ast+\bm{\xi}^\ast\cdot \bm{\nabla}^\ast f^\ast=-\frac{1}{\mathrm{Kn}}\left(f^\ast-{f^{eq}}^\ast\right),
 \end{equation}
 $$
-where the space, time, and microscopic velocity is omitted and where
+where the space, time, and microscopic velocity dependence is omitted and where
 $$
 \mathrm{Kn}=\frac{\tau\xi_0}{L},
 $$
diff --git a/sage/hermite.sage b/sage/hermite.sage
index 56b18ab..5a5af1d 100644
--- a/sage/hermite.sage
+++ b/sage/hermite.sage
@@ -1,16 +1,28 @@
 xi = var('xi')
 u = var('u')
 theta = var('theta')
+xi0 = var('xi0')
+n = var('n')
 assume(theta > 0)
 
-h0 = 1
-h1 = xi
-h2 = xi^2 - 1
-h3 = xi^3 - 3*xi
-h4 = xi^4 - 6*xi^2 + 3
-h5 = xi^5 - 10*xi^3 + 15*xi
+def hn(n):
+    return (-1)^n * xi0^n * exp(xi^2/(2 * xi0^2))*diff(exp(-xi^2/(2 * xi0^2)), xi, n)
 
-w = 1/sqrt(2*pi) * exp(-xi^2/2)
+print("Hermite polynomials")
+h0 = hn(0)
+print(h0)
+h1 = hn(1)
+print(h1)
+h2 = hn(2).full_simplify()
+print(h2)
+h3 = hn(3).full_simplify()
+print(h3)
+h4 = hn(4).full_simplify()
+print(h4)
+h5 = hn(5).full_simplify()
+print(h5)
+
+w = 1/sqrt(2*pi) * exp(-xi^2/(2*xi0^2))
 
 print("Coefficients")
 n0 = integrate(w * h0, xi, -infinity, infinity)
@@ -27,7 +39,7 @@ n5 = integrate(w * h5 * h5, xi, -infinity, infinity)
 print(n5.canonicalize_radical())
 
 
-bol = 1 / (sqrt(2 * pi * theta)) * exp(-(xi - u)^2/(2 * theta))
+bol = 1 / (sqrt(2 * pi * xi0^2)) * exp(-(xi - u)^2/(2 * xi0^2))
 res = integrate(bol, xi, -infinity, infinity)
 
 print("Moments of Boltzmann")
-- 
GitLab