#+PROPERTY: header-args :session gaussian_naive :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 utilities #+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/65039db7ebd200270fddddaf5c1fbd2d4737e189.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 = 1000 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( #+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/39028b895a317697621e8690d1e5a3188979c592.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/87e915528e9868083f4d250ead95c97db4f28988.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/aa7baf4e1909008bf3cdfac62aa71b85c32c343a.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/fb15573ff7e193fc4d2752e136d1a3970973ad0b.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/e7b5c7b9d95999a09ed1ae91d807e2d0c6194082.svg]] :END: #+begin_src jupyter-python plt.imshow(τ_r_alpha_star_G12.real) #+end_src #+RESULTS: :RESULTS: : [[file:./.ob-jupyter/b7c42a6558ae85d3bb0a0f0f66cfd9598c6d6874.svg]] :END: #+begin_src jupyter-python plt.imshow(τ_r_alpha_star_G12.imag) #+end_src #+RESULTS: :RESULTS: : [[file:./.ob-jupyter/92bcbc68cd7f5b94275bad7ba08e746677d9663a.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/e55e72e4f12b4de29c8e3b27b20751383dad8b91.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/160f021fbf643985f4243856170f6b2f03ea5a74.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/282d480c02843739748353addf8264ebadddf976.svg]] :END: #+end_src #+begin_src jupyter-python scipy.integrate.trapz(flow, t) #+end_src #+RESULTS: : 1.039473736501838 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/0e0bf4a15bdd12394587d63966b7a064782624d2.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/c19af3b4a1062b83834ca58fbde475f3c55fc771.svg]] :END: #+begin_src jupyter-python scipy.integrate.trapz(flow_alt, t) #+end_src #+RESULTS: : 1.1423457634247889 * 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/318a5c5ec950d216b888ed74ade4e3584f7ebb71.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: :RESULTS: \(x \mapsto \text{1.8908489286231014} + \text{15.192611898242761}\,x + \text{43.27633825204506}\,x^{2} + \text{49.327191029134}\,x^{3} + \text{28.124395600960767}\,x^{4} + \text{7.313444336316629}\,x^{5} + \text{1.0}\,x^{6}\) :END: 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: :RESULTS: \(x \mapsto \text{4.254410089401978} + \text{34.18337677104621}\,x + \text{99.26260999572447}\,x^{2} + \text{126.17879171379425}\,x^{3} + \text{106.55622835420678}\,x^{4} + \text{65.78244078584642}\,x^{5} + \text{30.374395600960767}\,x^{6} + \text{7.313444336316629}\,x^{7} + \text{1.0}\,x^{8}\) :END: 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: :RESULTS: \(x \mapsto \text{2.4651931085584264} + \text{22.18343348787121}\,x + \text{75.66089297865284}\,x^{2} + \text{112.84892451644214}\,x^{3} + \text{104.42196182727498}\,x^{4} + \text{65.80539139821734}\,x^{5} + \text{30.374395600960767}\,x^{6} + \text{7.313444336316629}\,x^{7} + \text{1.0}\,x^{8}\) :END: 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 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) return Ω * (R_l[None, :] * np.exp(t[:, None] * master_roots[None, :])).real.sum( axis=1 ) def G_11_ex(t): t = np.asarray(t) return ( R_l[None, :] ,* master_roots[None, :] ,* np.exp(t[:, None] * master_roots[None, :]) ).real.sum(axis=1) 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/c2fe807a440af4e45d937a16ab852638aa24bb3e.svg]] :END: #+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) def f_test(z): return ( Ω ** 2 + (Ω * (-G * ((z + γ) * s + δ * c) / (δ ** 2 + (γ + z) ** 2))).sum() + z ** 2 ) mpmath.findroot(f_test, Ω, solver="muller") #+end_src #+RESULTS: : mpc(real='-0.19348529127068425', imag='1.4906336228104743') #+begin_src jupyter-python plt.plot(τ, G_21_ex_alt(τ) - G_21_ex(τ)) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/8aff3059b259f06b7c7b98eedbb33c340037d0fb.svg]] :END: