master-thesis/python/energy_flow_proper/03_gaussian/analytic_calculation.org

• System Set-Up
• Exponential Fit
• Analytical Solution
  • Calculate the Magic Numbers
  • Construct the Polynomials
  • Calculate the Residuals
• Energy Flow

In this notebook, we compare the analytic solution with the numeric result. A nicer version with the code outsourced into a proper python module can be found here.

  import scipy
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy import integrate
  import sys
  import utilities
  import quadpy
  from hops.util import bcf
  import itertools

System Set-Up

  Ω = 1.5
  A = np.array([[0, Ω], [-Ω, 0]])
  η = 2
  ω_c = 1
  t_max = 25

  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

Exponential Fit

First we need an exponential fit for our BCF.

  W, G_raw = utilities.fit_α(α, 3, 80, 10_000)
  L, P = utilities.fit_α(α_0_dot, 5, 80, 10_000)

Looks quite good.

  τ = np.linspace(0, t_max, 1000)
  with utilities.hiro_style():
      fig, ax = utilities.plot_complex(τ, α(τ), label="exact")
      utilities.plot_complex(
          τ, utilities.α_apprx(τ, G_raw, W), ax=ax, label="fit", linestyle="--"
      )

  τ = np.linspace(0, -t_max, 1000)
  with utilities.hiro_style():
      fig, ax = utilities.plot_complex(τ, α(τ), label="exact")
      utilities.plot_complex(
          τ, utilities.α_apprx(-τ, G_raw.conj(), W.conj()), ax=ax, label="fit", linestyle="--"
      )

  fig, ax = utilities.plot_complex(τ, α_0_dot(τ), label="exact")
  utilities.plot_complex(
      τ, utilities.α_apprx(τ, P, L), ax=ax, label="fit", linestyle="--"
  )

hline

<AxesSubplot:>

Analytical Solution

Calculate the Magic Numbers

We begin with the $\varphi_n$ and $G_n$ from the original G.

  φ = -np.angle(G_raw)
  φ

array([-1.18710245,  0.64368323,  2.11154332])

  G = np.abs(G_raw)
  G

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$.

  γ, δ = W.real, W.imag

Now the \(s_n, c_n\).

  s, c = np.sin(φ), np.cos(φ)

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.

  roots_f = -(γ + δ * c/s)
  roots_f

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.

  roots_z = -W
  roots_z

array([-2.29230292-2.71252483j, -1.07996581-0.719112j  ,
       -0.28445344-0.09022829j])

Construct the Polynomials

  from numpy.polynomial import Polynomial

We later need \(f_0(z) = \prod_k (z-z_k) (z-z_k^{\ast})\).

  f_0 = utilities.poly_real(Polynomial.fromroots(np.concatenate((roots_z, roots_z.conj()))))
  f_0

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.

  p_1 = Polynomial([Ω**2, 0, 1]) * f_0
  p_1

The next ones are given through \(q_n=\Omega f_n(z) \prod_{k\neq n}(z-z_k) (z-z_k^{\ast})\)

  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)))
  ]

With this we construct our master polynomial \(p = p_1 + \sum_n q_n\).

  p = p_1 + sum(q)
  p

And find its roots.

  master_roots = p.roots()
  master_roots

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.

  np.unique(master_roots).size == master_roots.size

True

Very nice!

Calculate the Residuals

These are the prefactors for the diagonal.

  R_l = f_0(master_roots) / p.deriv()(master_roots)
  R_l

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\).

  def α_tilde(z):
      return (
          -G[None, :]
          ,* ((z[:, None] + γ[None, :]) * s[None, :] + δ[None, :] * c[None, :])
          / (δ[None, :] ** 2 + (γ[None, :] + z[:, None]) ** 2)
      ).sum(axis=1)

And these are for the most compliciated element.

  R_l_21 = (Ω + α_tilde(master_roots))* f_0(master_roots) / p.deriv()(master_roots)
  R_l_21

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\).

  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)
      )

  with utilities.hiro_style():
      plt.plot(τ, G_ex_new(τ).reshape(len(τ), 4))

Energy Flow

We adopt the terminology from my notes.

  A = G_11_ex
  B = G_12_ex
  Pk = P
  Bk = R_l * Ω
  Gk = G_raw
  Lk = L
  Ck = -master_roots # Mind the extra sign!
  Wk = W
  Ak = master_roots * R_l
  n = 1
  q_s_0 = 1 + 2 * n
  p_s_0 = 1 + 2 * n
  qp_0 = 1j

First, we calculate the coefficient matrices.

  Γ1 = (Bk[:, None] * Pk[None, :]) / (Lk[None, :] - Ck[:, None])
  Γ2 = (Bk[:, None] * Gk[None, :]) / (Ck[:, None] - Wk[None, :])
  Γ3 = (Bk[:, None] * Gk.conj()[None, :]) / (Ck[:, None] + Wk.conj()[None, :])
  ΓA = (Ak[:, None] * Pk[None, :]) / (Lk[None, :] - Ck[:, None])

The convolution free part is easy to calculate.

  def A_conv(t):
      t = np.asarray(t)
      result = np.zeros_like(t, dtype="complex128")
      for (n, m) in itertools.product(range(len(Ak)), range(len(Pk))):
          result += ΓA[n, m] * (np.exp(-Ck[n] * t) - np.exp(-Lk[m] * t))

      return result


  def B_conv(t):
      t = np.asarray(t)
      result = np.zeros_like(t, dtype="complex128")
      for (n, m) in itertools.product(range(len(Bk)), range(len(Pk))):
          result += Γ1[n, m] * (np.exp(-Ck[n] * t) - np.exp(-Lk[m] * t))

      return result


  def flow_conv_free(t):
      a, b = A(t), B(t)
      ac, bc = A_conv(t), B_conv(t)
      return (
          -1
          / 2
          ,* (
              q_s_0 * a * ac + p_s_0 * b * bc + qp_0 * a * bc + qp_0.conjugate() * b * ac
          ).imag
      )

The part containing the convolution with the BCF is a bit mode involved. The first version originates from my calculations and the second one arose from sagemath. Thank God! They do aggree now.

  def conv_part(t):
      t = np.asarray(t)
      result = np.zeros_like(t, dtype="float")
      for (m, k, n, l) in itertools.product(
          range(len(Bk)), range(len(Pk)), range(len(Bk)), range(len(Gk))
      ):
          g_1_2 = (
              Γ1[m, k]
              ,* Γ2[n, l]
              ,* (
                  (1 - np.exp(-(Ck[m] + Wk[l]) * t)) / (Ck[m] + Wk[l])
                  + (np.exp(-(Ck[m] + Ck[n]) * t) - 1) / (Ck[m] + Ck[n])
                  + (np.exp(-(Lk[k] + Wk[l]) * t) - 1) / (Lk[k] + Wk[l])
                  + (1 - np.exp(-(Lk[k] + Ck[n]) * t)) / (Lk[k] + Ck[n])
              )
          )

          g_1_3 = (
              Γ1[m, k]
              ,* Γ3[n, l]
              ,* (
                  (1 - np.exp(-(Ck[m] + Ck[n]) * t)) / (Ck[m] + Ck[n])
                  - (1 - np.exp(-(Lk[k] + Ck[n]) * t)) / (Lk[k] + Ck[n])
                  - (np.exp(-(Ck[n] + Wk[l].conj() * t)) - np.exp(-(Ck[m] + Ck[n]) * t))
                  / (Ck[m] - Wk[l].conj())
                  + (np.exp(-(Ck[n] + Wk[l].conj() * t)) - np.exp(-(Lk[k] + Ck[n]) * t))
                  / (Lk[k] - Wk[l].conj())
              )
          )

          result += -1 / 2 * (g_1_2.imag + g_1_3.imag)

      return result


  import monster_formula


  def conv_part_sage(t):
      t = np.asarray(t)
      result = np.zeros_like(t, dtype="float")
      for (n, k, l, m) in itertools.product(
          range(len(Bk)), range(len(Pk)), range(len(Gk)), range(len(Bk))
      ):

          result += (
              -1
              / 2
              ,* (
                  monster_formula.conv_part(
                      t,
                      Pk[k],
                      Lk[k],
                      Bk[n],
                      Ck[n],
                      Bk[m],
                      Ck[m],
                      Gk[l],
                      Wk[l],
                      Gk[l].conj(),
                      Wk[l].conj(),
                  )
              ).imag
          )

      return result

The contribution of the convolution part is pretty minor for zero temperature.

  τ = np.linspace(0, t_max, 5000)
  plt.plot(τ, conv_part(τ))

<matplotlib.lines.Line2D

at

0x7f85099ac190>

Now we can define the flow:

  def flow(t):
      return flow_conv_free(t) + conv_part(t)

And plot it against the numerically obtained flow.

  with open('res.npy', 'rb') as f:
      old_τ = np.load(f)
      old_flow_τ = np.load(f)
  with utilities.hiro_style():
      plt.plot(old_τ, flow(old_τ))
      plt.plot(old_τ, old_flow_τ)

They do aggree but I think the accuracy of the numerical integration is getting worse with time.

  scipy.integrate.quad(flow, 0, 20)

1.6004236143630652

1.564508138231918e-08

Integresting… maybe the fact that it does not converge to zero exactly has to do with this…