From e67cca834834f9bcb986560fc65e1c17602f6a53 Mon Sep 17 00:00:00 2001
From: Orestis <orestis.malaspinas@pm.me>
Date: Wed, 20 Mar 2024 22:36:28 +0100
Subject: [PATCH] started changing this to time-dependent stuff
---
sage/hermite.sage | 107 ++++++++++++++++++++++++++--------------------
1 file changed, 60 insertions(+), 47 deletions(-)
diff --git a/sage/hermite.sage b/sage/hermite.sage
index 5a5af1d..9db7805 100644
--- a/sage/hermite.sage
+++ b/sage/hermite.sage
@@ -1,12 +1,19 @@
+t = var('t')
+x = var('x')
+
xi = var('xi')
-u = var('u')
+u = function('u')(x, t)
theta = var('theta')
+rho = function('rho')(x, t)
+theta0 = var('theta0')
xi0 = var('xi0')
n = var('n')
+T = function('T')(x, t)
assume(theta > 0)
+assume(theta0 > 0)
def hn(n):
- return (-1)^n * xi0^n * exp(xi^2/(2 * xi0^2))*diff(exp(-xi^2/(2 * xi0^2)), xi, n)
+ return (-1)^n * theta0^(n/2) * exp(xi^2/(2 * theta0))*diff(exp(-xi^2/(2 * theta0)), xi, n)
print("Hermite polynomials")
h0 = hn(0)
@@ -22,7 +29,7 @@ print(h4)
h5 = hn(5).full_simplify()
print(h5)
-w = 1/sqrt(2*pi) * exp(-xi^2/(2*xi0^2))
+w = 1/sqrt(2*pi*theta0) * exp(-xi^2/(2*theta0))
print("Coefficients")
n0 = integrate(w * h0, xi, -infinity, infinity)
@@ -33,13 +40,13 @@ n2 = integrate(w * h2 * h2, xi, -infinity, infinity)
print(n2.canonicalize_radical())
n3 = integrate(w * h3 * h3, xi, -infinity, infinity)
print(n3.canonicalize_radical())
-n4 = integrate(w * h4 * h4, xi, -infinity, infinity)
-print(n4.canonicalize_radical())
-n5 = integrate(w * h5 * h5, xi, -infinity, infinity)
-print(n5.canonicalize_radical())
+# n4 = integrate(w * h4 * h4, xi, -infinity, infinity)
+# print(n4.canonicalize_radical())
+# n5 = integrate(w * h5 * h5, xi, -infinity, infinity)
+# print(n5.canonicalize_radical())
-bol = 1 / (sqrt(2 * pi * xi0^2)) * exp(-(xi - u)^2/(2 * xi0^2))
+bol = rho / (sqrt(2 * pi * theta)) * exp(-(xi - u)^2/(2 * theta))
res = integrate(bol, xi, -infinity, infinity)
print("Moments of Boltzmann")
@@ -51,47 +58,53 @@ mom2 = integrate(bol * h2, xi, -infinity, infinity)
print(mom2.canonicalize_radical())
mom3 = integrate(bol * h3, xi, -infinity, infinity)
print(mom3.canonicalize_radical())
-mom4 = integrate(bol * h4, xi, -infinity, infinity)
-print(mom4.canonicalize_radical())
-mom5 = integrate(bol * h5, xi, -infinity, infinity)
-print(mom5.canonicalize_radical())
+# mom4 = integrate(bol * h4, xi, -infinity, infinity)
+# print(mom4.canonicalize_radical())
+# mom5 = integrate(bol * h5, xi, -infinity, infinity)
+# print(mom5.canonicalize_radical())
-print("Moments of Polynomial")
-feq = w * (mom0 + h1 * mom1/n1 + h2 * mom2/n2 + h3 * mom3/n3 + h4 * mom4/n4 + h5 * mom5/n5)
-p_mom0 = integrate(feq * 1, xi, -infinity, infinity)
-print(p_mom0.canonicalize_radical())
-p_mom1 = integrate(feq * h1, xi, -infinity, infinity)
-print(p_mom1.canonicalize_radical())
-p_mom2 = integrate(feq * h2, xi, -infinity, infinity)
-print(p_mom2.canonicalize_radical())
-p_mom3 = integrate(feq * h3, xi, -infinity, infinity)
-print(p_mom3.canonicalize_radical())
-p_mom4 = integrate(feq * h4, xi, -infinity, infinity)
-print(p_mom4.canonicalize_radical())
-p_mom5 = integrate(feq * h5, xi, -infinity, infinity)
-print(p_mom5.canonicalize_radical())
+print(diff((integrate(bol * h0 * xi, xi, -infinity, infinity).canonicalize_radical() / sqrt(theta0)), t))
+print((integrate(bol * h1 * xi, xi, -infinity, infinity).canonicalize_radical() / sqrt(theta0)).subs(theta=T))
+print((integrate(bol * h2 * xi, xi, -infinity, infinity).canonicalize_radical() / sqrt(theta0)).subs(theta=T))
-w_m = var('w_m')
-w_0 = var('w_0')
-w_p = var('w_p')
-xi_m = var('xi_m')
-xi_0 = var('xi_0')
-xi_p = var('xi_p')
-wi = vector([-w_p, w_0, w_p])
-xi_i = vector([-xi_p, xi_0, xi_p])
-
-h0_i = 1
-h1_i = xi_i
-h2_i = vector([xi_i[i]^2 - 1 for i in range(len(xi_i))])
-h3_i = vector([xi_i[i]^3 - 3 * xi_i[i] for i in range(len(xi_i))])
-h4_i = vector([xi_i[i]^4 - 6 * xi_i[i]^2 + 3 for i in range(len(xi_i))])
-h5_i = vector([xi_i[i]^5 - 10 * xi_i[i]^3 + 15 * xi_i[i] for i in range(len(xi_i))])
-sum(wi) == n0
-wi.dot_product(h1_i) == n1
-wi.dot_product(h2_i) == n2
-wi.dot_product(h3_i) == n3
-wi.dot_product(h4_i) == n4
-wi.dot_product(h5_i) == n5
+# print("Moments of Polynomial")
+# feq = w * (mom0 + h1 * mom1/n1 + h2 * mom2/n2 + h3 * mom3/n3 + h4 * mom4/n4 + h5 * mom5/n5)
+# p_mom0 = integrate(feq * 1, xi, -infinity, infinity)
+# print(p_mom0.canonicalize_radical())
+# p_mom1 = integrate(feq * h1, xi, -infinity, infinity)
+# print(p_mom1.canonicalize_radical())
+# p_mom2 = integrate(feq * h2, xi, -infinity, infinity)
+# print(p_mom2.canonicalize_radical())
+# p_mom3 = integrate(feq * h3, xi, -infinity, infinity)
+# print(p_mom3.canonicalize_radical())
+# p_mom4 = integrate(feq * h4, xi, -infinity, infinity)
+# print(p_mom4.canonicalize_radical())
+# p_mom5 = integrate(feq * h5, xi, -infinity, infinity)
+# print(p_mom5.canonicalize_radical())
+# w_m = var('w_m')
+# w_0 = var('w_0')
+# w_p = var('w_p')
+# xi_m = var('xi_m')
+# xi_0 = var('xi_0')
+# xi_p = var('xi_p')
+# wi = vector([-w_p, w_0, w_p])
+# xi_i = vector([-xi_p, xi_0, xi_p])
+#
+# h0_i = 1
+# h1_i = xi_i
+# h2_i = vector([xi_i[i]^2 - 1 for i in range(len(xi_i))])
+# h3_i = vector([xi_i[i]^3 - 3 * xi_i[i] for i in range(len(xi_i))])
+# h4_i = vector([xi_i[i]^4 - 6 * xi_i[i]^2 + 3 for i in range(len(xi_i))])
+# h5_i = vector([xi_i[i]^5 - 10 * xi_i[i]^3 + 15 * xi_i[i] for i in range(len(xi_i))])
+#
+# sum(wi) == n0
+# wi.dot_product(h1_i) == n1
+# wi.dot_product(h2_i) == n2
+# wi.dot_product(h3_i) == n3
+# wi.dot_product(h4_i) == n4
+# wi.dot_product(h5_i) == n5
+#
+#
--
GitLab