#+PROPERTY: header-args :session gaussian :kernel python :pandoc yes :async yes #+begin_src jupyter-python :results none import scipy import numpy as np import matplotlib.pyplot as plt from scipy import integrate import sys sys.path.append("/home/hiro/python/energy_flow_proper/hops") sys.path.append("/home/hiro/python/energy_flow_proper") import utilities import quadpy #+end_src * System Parameters #+begin_src jupyter-python :results none Ω = 1.5 A = np.array([[0, Ω], [-Ω, 0]]) η = 2 ω_c = 1 t_max = 16 def α(t): return η / np.pi * (ω_c / (1 + 1j * ω_c * t)) ** 2 def α_0_dot(t): return 2 * η / (1j * np.pi) * (ω_c / (1 + 1j * ω_c * t)) ** 3 def K(τ): return np.array([[0, 0], [α(τ).imag, 0]]) def L(t, s): return np.exp(-A * t) @ K(t - s) @ np.exp(A * s) #+end_src * Definition of the RHS #+begin_src jupyter-python :results none def H_dot(t, _): return -integrate.trapezoid( [(L(t, ts[s]) @ np.array(H_to_date[s])).flatten() for s in range(len(ts))], ts, axis=0, ) def K_dot(t, curr): return (A @ curr.reshape(2, 2)).flatten() - integrate.simpson( [(K(t - ts[s]) @ np.array(H_to_date[s])).flatten() for s in range(len(ts))] + [(K(0) @ curr.reshape(2, 2)).flatten()], ts + [t], axis=0, ) #+end_src * Scipy Hack Hacking together an integro-differential solver. #+begin_src jupyter-python :results none ts = [0] H_to_date = [np.eye(2)] r = integrate.ode(K_dot).set_integrator('vode') r.set_initial_value(np.eye(2, dtype="float").flatten(), 0) t1 = t_max dt = .001 while r.successful() and r.t < t1: sol = r.integrate(r.t + dt) H_to_date.append(sol.reshape(2,2)) ts.append(r.t + dt) #+end_src Dividing out the unitary dynamics. #+begin_src jupyter-python #plt.plot(ts, (np.array([scipy.linalg.expm(-A * np.array(t)) for t in ts]) @ np.array(H_to_date))[:,1,0]) proper = np.array(H_to_date) #+end_src #+RESULTS: \(\langle q\rangle\) #+begin_src jupyter-python plt.plot(ts, proper @ np.array([1, 0])) #+end_src #+RESULTS: :RESULTS: | | | [[file:./.ob-jupyter/699a32087aa1ac7105aef99a1b81469bbeb8d28a.svg]] :END: This looks like a slightly dampened HO. * IDESolver With a proper solver it works like a charm. See [[https://idesolver.readthedocs.io/en/latest/manual.html][the docs]]. #+begin_src jupyter-python from idesolver import IDESolver steps = 500 t = np.linspace(0, t_max, steps) solver = IDESolver( x = t, y_0 = np.array([1, 0, 0, 1]), c = lambda x, y: (A @ y.reshape(2, 2)).flatten(), d = lambda x: 1, k = lambda x, s: -K(x-s)[1, 0], f = lambda y: np.array([0, 0, y[0], y[1]]), lower_bound = lambda x: 0, upper_bound = lambda x: x, ) sol = solver.solve() #+end_src #+RESULTS: #+begin_example /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 0 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 1 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 2 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 3 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 4 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 5 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 6 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 7 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 8 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 9 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 10 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 11 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 12 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 13 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 14 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 15 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 16 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 17 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 18 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 19 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 20 warnings.warn( /nix/store/j2851gf1ssx8b9qz14b8n44f1cq5dhi9-python3-3.9.4-env/lib/python3.9/site-packages/idesolver/idesolver.py:292: IDEConvergenceWarning: Error increased on iteration 24 warnings.warn( #+end_example Reshape \(G\) into a time series of matrices. #+begin_src jupyter-python G = np.einsum("ijk->kij", sol.reshape((2,2,steps))) #+end_src #+RESULTS: And plot the time development of the mean values of \(p\) and \(p\). #+begin_src jupyter-python plt.plot(t, G @ np.array([1, 0])) #+end_src #+RESULTS: :RESULTS: | | | [[file:./.ob-jupyter/fad23f518c037f04e6108539af9a0652901fe1ff.svg]] :END: Looks OK. #+begin_src jupyter-python plt.plot( t, (np.array([scipy.linalg.expm(-A * np.array(τ)) for τ in t]) @ G).reshape(steps, 4), ) #+end_src #+RESULTS: :RESULTS: | | | | | [[file:./.ob-jupyter/e731bf5509fa205c540e9062c3a15c878b8a4a21.svg]] :END: We see that the derivative at the begining is really 0 if we divide out the unitary dynamics. The scipy integration arrives at pretty much the same picture :). #+begin_src jupyter-python plt.plot( ts, (np.array([scipy.linalg.expm(-A * np.array(t)) for t in ts]) @ np.array(H_to_date)).reshape(len(ts), 4), label="scipy" ) plt.plot( t, (np.array([scipy.linalg.expm(-A * np.array(τ)) for τ in t]) @ G).reshape(steps, 4), linestyle="--", label="idesolve" ) plt.legend() plt.xlabel('t') plt.ylabel('$e^{-At}G_{ij}$') #+end_src #+RESULTS: :RESULTS: : Text(0, 0.5, '$e^{-At}G_{ij}$') [[file:./.ob-jupyter/18f4e12520e5e460dfc0d1a5b4ba87d7d087be39.svg]] :END: The initial slip is removed in idesolve due to the iteration. You have to look at the off diagonal curves to see that. #+begin_src jupyter-python plt.plot(t, np.linalg.det(G)) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/0d736728dd6cac1d188363bc857527c2c70fcc4e.svg]] :END: We see that the "volume" bleeds out of the system. * Interpolation #+begin_src jupyter-python import fcSpline #+end_src #+RESULTS: #+begin_src jupyter-python class FCSWrap: def __init__(self, t_min, t_max, G): self._G = G self._t_min = t_min self._t_max = t_max self._splines = np.array( list(map(lambda g: fcSpline.FCS(self._t_min, self._t_max, g), G)) ) def __call__(self, t): t = np.asarray(t).astype('float64') res = np.array(list(map(lambda g: g(t), self._splines))) if t.size == 1: return res.reshape((2,2)) return res.reshape((2,2,t.size)).swapaxes(0,2).swapaxes(1,2) def __getitem__(self, key): return self._splines.reshape(2,2)[key] #+end_src #+RESULTS: #+begin_src jupyter-python G_inter = FCSWrap(t.min(), t.max(), sol) #+end_src #+RESULTS: * Calculate the Convolutions We use quadpy for complex integration. #+begin_src jupyter-python import quadpy #+end_src #+RESULTS: We will need \(G_{12}\) often. #+begin_src jupyter-python G12 = G_inter[0, 1] #+end_src #+RESULTS: Calculate \(G_{12}\ast \dot{\alpha_0}\). #+begin_src jupyter-python :results none G12_star_α_0_dot = np.array( [ quadpy.quad(lambda τ: G12(τ) * α_0_dot(t_now - τ), 0, t_now)[0] for t_now in t ] ) G12_star_α_0_dot = fcSpline.FCS(t.min(), t.max(), G12_star_α_0_dot) #+end_src Calculate \(\tau_r\alpha\ast G_{12}\). #+begin_src jupyter-python :results none τ_r_alpha_star_G12 = np.array( [ [ quadpy.quad(lambda τ: G_inter[0, 1](t_now - τ) * α_0_dot(τ - r), 0, t_now)[ 0 ] for r in t ] for t_now in t ] ) τ_r_alpha_star_G12_real = scipy.interpolate.interp2d(t, t, τ_r_alpha_star_G12.real) τ_r_alpha_star_G12_imag = scipy.interpolate.interp2d(t, t, τ_r_alpha_star_G12.imag) def τ_r_alpha_star_G12_inter(t, r): return τ_r_alpha_star_G12_real(t, r) + 1j * τ_r_alpha_star_G12_imag(t, r) #+end_src #+begin_src jupyter-python plt.plot(t, G12_star_α_0_dot(t).imag) plt.plot(t, G12_star_α_0_dot(t).real) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/71d347c687e062a6feb23ab0f42f6b4f328155b6.svg]] :END: #+begin_src jupyter-python plt.imshow(τ_r_alpha_star_G12.real) #+end_src #+RESULTS: :RESULTS: : [[file:./.ob-jupyter/0dcad0dda95f3d0c88222bfc06e3ab453b898aa5.svg]] :END: #+begin_src jupyter-python plt.imshow(τ_r_alpha_star_G12.imag) #+end_src #+RESULTS: :RESULTS: : [[file:./.ob-jupyter/1c3acd3be3596cb399b651d18fe9b3777dcdf185.svg]] :END: And the final integral. #+begin_src jupyter-python flow_half = np.array( [ quadpy.quad( lambda r: G12_star_α_0_dot(t_now - r) * τ_r_alpha_star_G12_inter(t_now, r)[0], 0, t_now, )[0] for t_now in t ] ) flow_half = -1/2 * flow_half.imag #+end_src #+RESULTS: #+begin_src jupyter-python plt.plot(t, flow_half) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/3fd6e2e90774a32ecd7037d93121032e347c999f.svg]] :END: And now the \(G\) part of the flow. #+begin_src jupyter-python n = 1 flow_G_half = np.array( [ quadpy.quad( lambda s: ( ( G_inter[0, 0](t_now) * G_inter[0, 1](s) - G_inter[0, 0](s) * G_inter[0, 1](t_now) ) ,* α_0_dot(t_now - s).real + (2 * n + 1) ,* ( G_inter[0, 0](t_now) * G_inter[0, 0](s) + G_inter[0, 1](t_now) * G_inter[0, 1](s) ) ,* α_0_dot(t_now - s).imag ), 0, t_now, )[0] for t_now in t ] ) flow_G_half = -1 / 2 * flow_G_half #+end_src #+RESULTS: #+begin_src jupyter-python plt.plot(t, flow_G_half) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/df69413ad6fb4f37ec6ad22639acc3cebdf7ee0b.svg]] :END: #+begin_src jupyter-python flow = flow_G_half + flow_half #+end_src #+RESULTS: #+begin_src jupyter-python plt.plot(t, flow) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/88b32148be473f0d22726d145bbacf7a6acdf794.svg]] :END: #+end_src #+begin_src jupyter-python scipy.integrate.trapz(flow, t) #+end_src #+RESULTS: : 1.5215449446799811 There remains the zero point energy. #+begin_src jupyter-python def integrand(r, s, l, t_now): return -1 / 2 * (G12(t_now - l) * G12(s - r) * α(l - r) * α_0_dot(t_now - s)).imag flow_half_alt = np.array( [ scipy.integrate.nquad( integrand, [lambda s, _, __: [0, s], [0, t_now], [0, t_now]], args=[t_now] ) for t_now in np.linspace(0, t.max(), 100) ] ) #+end_src #+RESULTS: #+begin_src jupyter-python flow_half_alt_int = fcSpline.FCS(0, t.max(), flow_half_alt[:, 0]) #+end_src #+RESULTS: #+begin_src jupyter-python plt.plot(t, flow_half_alt_int(t)) plt.plot(t, flow_half) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/1d9aaf9144c9a0f5e43b0c97e9b8d1966188f3a8.svg]] :END: #+begin_src jupyter-python flow_alt = flow_G_half + flow_half_alt_int(t) #+end_src #+RESULTS: #+begin_src jupyter-python plt.plot(t, flow_alt) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/8a927282d95f0cc9b714cba88fa99ed00ecf5f4f.svg]] :END: #+begin_src jupyter-python scipy.integrate.trapz(flow_alt, t) #+end_src #+RESULTS: : 1.5885826683889948 * Analytical Solution ** Exponential Fit First we need an exponential fit for our BCF. #+begin_src jupyter-python W, G_raw = utilities.fit_α(α, 3, 80, 10_000) τ = np.linspace(0, t_max, 1000) #+end_src #+RESULTS: Looks quite good. #+begin_src jupyter-python fig, ax = utilities.plot_complex(τ, α(τ), label="exact") utilities.plot_complex( τ, utilities.α_apprx(τ, G_raw, W), ax=ax, label="fit", linestyle="--" ) #+end_src #+RESULTS: :RESULTS: | hline | | [[file:./.ob-jupyter/afa5b41a1ab07e555bddc54be560b06c11bf3906.svg]] :END: ** Calculate the Magic Numbers We begin with the $\varphi_n$ and $G_n$ from the original ~G~. #+begin_src jupyter-python φ = -np.angle(G_raw) φ #+end_src #+RESULTS: : array([-1.18710245, 0.64368323, 2.11154332]) #+begin_src jupyter-python G = np.abs(G_raw) G #+end_src #+RESULTS: : array([0.51238635, 0.62180167, 0.10107935]) Now we calculate the real and imag parts of the ~W~ parameters and call them $\gamma_n$ and $\delta_n$. #+begin_src jupyter-python γ, δ = W.real, W.imag #+end_src #+RESULTS: Now the \(s_n, c_n\). #+begin_src jupyter-python s, c = np.sin(φ), np.cos(φ) #+end_src #+RESULTS: Now we calculate the roots of $f_n(z)=-G_n ((z+\gamma_n) s_n + \delta_n c_n)$. Normally we should be vary of one of the \(\deltas\) being zero, but this is not the case. #+begin_src jupyter-python roots_f = -(γ + δ * c/s) roots_f #+end_src #+RESULTS: : array([-1.19725058, -2.0384181 , -0.23027627]) Now the \(z_k\) the roots of \(\delta_k^2 + (\gamma_k + z)^2\). *We don't include the conjugates.* #+begin_src jupyter-python roots_z = -W roots_z #+end_src #+RESULTS: : array([-2.29230292-2.71252483j, -1.07996581-0.719112j , : -0.28445344-0.09022829j]) ** Construct the Polynomials #+begin_src jupyter-python from numpy.polynomial import Polynomial #+end_src #+RESULTS: We later need \(f_0(z) = \prod_k (z-z_k) (z-z_k^{\ast})\). #+begin_src jupyter-python f_0 = utilities.poly_real(Polynomial.fromroots(np.concatenate((roots_z, roots_z.conj())))) f_0 #+end_src #+RESULTS: Another polynomial is simply \(p_1 = (z^2 + \Omega^2)\prod_k (z-z_k) (z-z_k^{\ast})\) and we can construct it from its roots. #+begin_src jupyter-python p_1 = Polynomial([Ω**2, 0, 1]) * f_0 p_1 #+end_src #+RESULTS: The next ones are given through \(q_n=\Omega f_n(z) \prod_{k\neq n}(z-z_k) (z-z_k^{\ast})\) #+begin_src jupyter-python q = [ -G_c ,* Ω * s_c ,* utilities.poly_real(Polynomial.fromroots( np.concatenate( ( [root_f], utilities.except_element(roots_z, i), utilities.except_element(roots_z, i).conj(), ) ) )) for root_f, G_c, γ_c, δ_c, s_c, c_c, i in zip(roots_f, G, γ, δ, s, c, range(len(c))) ] #+end_src #+RESULTS: With this we construct our master polynomial \(p = p_1 + \sum_n q_n\). #+begin_src jupyter-python p = p_1 + sum(q) p #+end_src #+RESULTS: And find its roots. #+begin_src jupyter-python master_roots = p.roots() master_roots #+end_src #+RESULTS: : array([-2.28139877-2.68284887j, -2.28139877+2.68284887j, : -0.93979297-0.62025112j, -0.93979297+0.62025112j, : -0.24204514-0.10390896j, -0.24204514+0.10390896j, : -0.19348529-1.49063362j, -0.19348529+1.49063362j]) Let's see if they're all unique. This should make things easier. #+begin_src jupyter-python np.unique(master_roots).size == master_roots.size #+end_src #+RESULTS: : True Very nice! ** Calculate the Residuals These are the prefactors for the diagonal. #+begin_src jupyter-python R_l = f_0(master_roots) / p.deriv()(master_roots) R_l #+end_src #+RESULTS: : array([ 0.00251589-0.00080785j, 0.00251589+0.00080785j, : 0.06323421+0.02800137j, 0.06323421-0.02800137j, : 0.02732669-0.00211665j, 0.02732669+0.00211665j, : -0.09307679+0.36145133j, -0.09307679-0.36145133j]) And the laplace transform of \(\alpha\). #+begin_src jupyter-python def α_tilde(z): return ( -G[None, :] ,* ((z[:, None] + γ[None, :]) * s[None, :] + δ[None, :] * c[None, :]) / (δ[None, :] ** 2 + (γ[None, :] + z[:, None]) ** 2) ).sum(axis=1) #+end_src #+RESULTS: And these are for the most compliciated element. #+begin_src jupyter-python R_l_21 = (Ω + α_tilde(master_roots))* f_0(master_roots) / p.deriv()(master_roots) R_l_21 #+end_src #+RESULTS: : array([-0.00325014-0.02160514j, -0.00325014+0.02160514j, : 0.00074814-0.05845205j, 0.00074814+0.05845205j, : -0.00094159-0.00084894j, -0.00094159+0.00084894j, : 0.00344359+0.56219924j, 0.00344359-0.56219924j]) Now we can calculate \(G\). #+begin_src jupyter-python def G_12_ex(t): t = np.asarray(t) t_shape = t.shape if len(t.shape) == 0: t = t.reshape((1,)) return Ω * (R_l[None, :] * np.exp(t[:, None] * master_roots[None, :])).real.sum( axis=1 ).reshape(t_shape) def G_11_ex(t): t = np.asarray(t) t_shape = t.shape if len(t.shape) == 0: t = np.array([t]) return ( R_l[None, :] ,* master_roots[None, :] ,* np.exp(t[:, None] * master_roots[None, :]) ).real.sum(axis=1).reshape(t_shape) def G_21_ex(t): t = np.asarray(t) return -(R_l_21[None, :] * np.exp(t[:, None] * master_roots[None, :])).real.sum( axis=1 ) def G_21_ex_alt(t): t = np.asarray(t) return ( R_l[None, :] ,* master_roots[None, :] ** 2 ,* np.exp(t[:, None] * master_roots[None, :]) ).real.sum(axis=1) / Ω def G_ex(t): t = np.asarray(t) if t.size == 1: t = np.array([t]) diag = G_11_ex(t) return ( np.array([[diag, G_12_ex(t)], [G_21_ex(t), diag]]).swapaxes(0, 2).swapaxes(1, 2) ) def G_ex_new(t): t = np.asarray(t) if t.size == 1: t = np.array([t]) g_12 = R_l[None, :] * np.exp(t[:, None] * master_roots[None, :]) diag = master_roots[None, :] * g_12 g_21 = master_roots[None, :] * diag / Ω return ( np.array([[diag, g_12 * Ω], [g_21, diag]]) .real.sum(axis=3) .swapaxes(0, 2) .swapaxes(1, 2) ) #+end_src #+RESULTS: #+begin_src jupyter-python #plt.plot(τ, G_ex(τ).reshape(len(τ), 4)) plt.plot(τ, G_ex_new(τ).reshape(len(τ), 4)) #plt.plot(ts, proper.reshape(len(ts), 4), linestyle="--") #+end_src #+RESULTS: :RESULTS: | | | | | [[file:./.ob-jupyter/caa8f4b42f6de38f8733ac5169a063de1d5c5d53.svg]] :END: #+begin_src jupyter-python plt.plot(τ, G_21_ex_alt(τ) - G_21_ex(τ)) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/5cf02217f725f3e631bb65db0df160a7584ebf65.svg]] :END: * Correlation Function We use the terminology form my notes. #+begin_src jupyter-python def B_s(t): return Ω * (R_l * np.exp(t * master_roots)).real.sum() def A_s(t): return ( R_l ,,* master_roots ,,* np.exp(t * master_roots) ).real.sum() A = G_11_ex B = G_12_ex n = 1 q_s_0 = 1 + 2 * n p_s_0 = 1 + 2 * n qp_0 = 1j #+end_src #+RESULTS: We'd like to calculate \(\langle q(t), q(s)\rangle\). The convolution free part is easy. #+begin_src jupyter-python def q_corr_conv_free(t, s): return ( q_s_0 * A(t) * A(s) + p_s_0 * B(t) * B(s) + qp_0 * A(t) * B(s) + qp_0.conjugate() * B(t) * A(s) ) #+end_src #+RESULTS: Now the convolution part (we don't really need this actually). #+begin_src jupyter-python def Λ_real(t, s): def integrand(l, r): return (B(t - l) * B(s - r) * α(l - r)).real return scipy.integrate.dblquad( integrand, 0, s, lambda _: 0, lambda _: t )[0] def Λ_real_mag(t, s): def integrand(l, r): return (B(t - l) * B(s - r) * α(l - r)).imag return scipy.integrate.dblquad( integrand, 0, s, lambda _: 0, lambda _: t )[0] #+end_src #+RESULTS: With this we can define the correlation function. #+begin_src jupyter-python def q_ts(t, s): return q_corr_conv_free(t, s) + Λ(t, s) #+end_src #+RESULTS: For the energy flow we have to calculate the convolution with \(\dot{\alpha_0}\). We can reduce the triple integral to a double integral by assuming the zero temp BCF to be a sum of exponentials. * Energy Flow This is "just" another convolution. We have to integrate over a wedge. #+begin_src jupyter-python from functools import cache, wraps @cache def flow(t): conv_free = quadpy.quad( lambda s: (q_corr_conv_free(t, s) * α_0_dot(t - s)).imag, 0, t )[0] def conv_integrand(l, r, s): return (B(t - l) * B(s - r) * α(l - r) * α_0_dot(t - s)).imag # def conv_integrand_quadpy(x): # s, r, l = x # if (s 2 from multiprocessing import Pool 3 with Pool(4) as p: ----> 4 flow_τ = np.array(p.map(flow, τ)) ~/.conda/envs/hops/lib/python3.9/multiprocessing/pool.py in map(self, func, iterable, chunksize) 362 in a list that is returned. 363 ''' --> 364 return self._map_async(func, iterable, mapstar, chunksize).get() 365 366 def starmap(self, func, iterable, chunksize=None): ~/.conda/envs/hops/lib/python3.9/multiprocessing/pool.py in get(self, timeout) 763 764 def get(self, timeout=None): --> 765 self.wait(timeout) 766 if not self.ready(): 767 raise TimeoutError ~/.conda/envs/hops/lib/python3.9/multiprocessing/pool.py in wait(self, timeout) 760 761 def wait(self, timeout=None): --> 762 self._event.wait(timeout) 763 764 def get(self, timeout=None): ~/.conda/envs/hops/lib/python3.9/threading.py in wait(self, timeout) 572 signaled = self._flag 573 if not signaled: --> 574 signaled = self._cond.wait(timeout) 575 return signaled 576 ~/.conda/envs/hops/lib/python3.9/threading.py in wait(self, timeout) 310 try: # restore state no matter what (e.g., KeyboardInterrupt) 311 if timeout is None: --> 312 waiter.acquire() 313 gotit = True 314 else: KeyboardInterrupt: #+end_example :END: #+begin_src jupyter-python plt.plot(τ, flow_τ) #+end_src #+RESULTS: #+begin_src jupyter-python scipy.integrate.trapz(flow_τ, τ) #+end_src #+RESULTS: #+begin_src jupyter-python with open("res.npy", 'wb') as f: np.save(f, τ) np.save(f, flow_τ) #+end_src #+begin_src jupyter-python with open("res.npy", 'rb') as f: τ = np.load(f) flow_τ = np.load(f) #+end_src #+RESULTS: #+begin_src jupyter-python with utilities.hiro_style(): plt.plot(τ, flow_τ) #+end_src #+RESULTS: [[file:./.ob-jupyter/f8925d0fb5e2dcf292b99f7784b208b718bc92cd.svg]] * Scratch #+begin_src jupyter-python from hops.util import bcf α_nonzero = bcf.OhmicBCF_nonZeroTemp(1, 1, 5, 1/2) W, G_raw = utilities.fit_α(α_nonzero, 6, 150, 10_000) τ = np.linspace(0, t_max, 1000) #+end_src #+RESULTS: Looks quite good. #+begin_src jupyter-python fig, ax = utilities.plot_complex(τ, α_nonzero(τ), label="exact") utilities.plot_complex( τ, utilities.α_apprx(τ, G_raw, W), ax=ax, label="fit", linestyle="--" ) #+end_src #+RESULTS: :RESULTS: | hline | | [[file:./.ob-jupyter/f0e3aad0ce345f9d66f2435a7f1dd6f85afa128b.svg]] :END: