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

  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)

  τ = np.linspace(0, t_max, 1000)

Looks quite good.

  fig, ax = utilities.plot_complex(τ, α(τ), label="exact")
  utilities.plot_complex(
      τ, utilities.α_apprx(τ, G_raw, W), ax=ax, label="fit", linestyle="--"
  )

hline

<AxesSubplot:>

hline

<AxesSubplot:>

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

  plt.plot(τ, G_ex_new(τ).reshape(len(τ), 4))

<matplotlib.lines.Line2D

at

0x7fef38be9460>

<matplotlib.lines.Line2D

at

0x7fef38be9490>

<matplotlib.lines.Line2D

at

0x7fef38be95b0>

<matplotlib.lines.Line2D

at

0x7fef38be96d0>

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 the coefficient matrices.

  # def Γ1(n, k):
  #     return (Bk[n] * Pk[k]) / (Lk[k] - Ck[n])]

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

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


  @np.vectorize
  def flow_conv_free_old(t):
      conv_free = quadpy.quad(
          lambda s: (q_corr_conv_free(t, s) * α_0_dot(t - s)).imag, 0, t
      )[0]

      return -1 / 2 * conv_free


  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
      )

Then the auxiliary integrals.

  import numba

  def conv_part(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(Bk)), range(len(Gk))
      ):
          g_1_2 = (
              Γ1[n, k]
              ,* Γ2[l, m]
              ,* (
                  (1 - np.exp(-(Ck[n] + Wk[m]) * t)) / (Ck[n] + Wk[m])
                  + (np.exp(-(Ck[n] + Ck[l]) * t) - 1) / (Ck[n] + Ck[l])
                  + (np.exp(-(Lk[k] + Wk[m]) * t) - 1) / (Lk[k] + Wk[m])
                  + (1 - np.exp(-(Lk[k] + Ck[l]) * t)) / (Lk[k] + Ck[l])
              )
          )

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

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

      return result

Seems to be correct :P or at least plausible.

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

<matplotlib.lines.Line2D

at

0x7fef12c0b190>

Now we can define the flow:

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

  plt.plot(τ, flow(τ))

<matplotlib.lines.Line2D

at

0x7fef12f25f70>

  scipy.integrate.trapz(flow(τ), τ)

1.7594404753272874