added sage file

sage/hermite.sage

+xi = var('xi')
+u = var('u')
+theta = var('theta')
+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
+
+w = 1/sqrt(2*pi) * exp(-xi^2/2)
+
+print("Coefficients")
+n0 = integrate(w * h0, xi, -infinity, infinity)
+print(n0.canonicalize_radical())
+n1 = integrate(w * h1 * h1, xi, -infinity, infinity)
+print(n1.canonicalize_radical())
+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())
+
+
+bol = 1 / (sqrt(2 * pi * theta)) * exp(-(xi - u)^2/(2 * theta))
+res = integrate(bol, xi, -infinity, infinity)
+
+print("Moments of Boltzmann")
+mom0 = integrate(bol * 1, xi, -infinity, infinity)
+print(mom0.canonicalize_radical())
+mom1 = integrate(bol * h1, xi, -infinity, infinity)
+print(mom1.canonicalize_radical())
+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())
+
+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
+
+