mirror of
https://github.com/vale981/master-thesis
synced 2025-03-09 04:36:39 -04:00
107 lines
2.5 KiB
Org Mode
107 lines
2.5 KiB
Org Mode
|
#+PROPERTY: header-args :session gaussian_farts :kernel python :pandoc yes :async yes
|
|||
|
|
|
|||
|
|
#+begin_src jupyter-python
|
|||
|
|
import scipy
|
|||
|
|
import numpy as np
|
|||
|
|
import matplotlib.pyplot as plt
|
|||
|
|
|
|||
|
|
|
|||
|
|
from scipy.special import expi
|
|||
|
|
#+end_src
|
|||
|
|
|
|||
|
|
#+RESULTS:
|
|||
|
|
|
|||
|
|
* Inverse Laplace Transform of the Propagator
|
|||
|
|
We need the laplace transform of our spectral density.
|
|||
|
|
#+begin_src jupyter-python
|
|||
|
|
η = 1
|
|||
|
|
ω_c = 1
|
|||
|
|
Ω = 1
|
|||
|
|
|
|||
|
|
def α(t):
|
|||
|
|
return η / np.pi * (ω_c / (1 + 1j * ω_c * t)) ** 2
|
|||
|
|
|
|||
|
|
|
|||
|
|
def α_l(z):
|
|||
|
|
return (
|
|||
|
|
(1j * ω_c * η)
|
|||
|
|
/ np.pi
|
|||
|
|
,* (1 + 1j * z / ω_c * np.exp(-1j * z / ω_c) * expi(1j * z / ω_c))
|
|||
|
|
)
|
|||
|
|
def α_l_num(z):
|
|||
|
|
return quadpy.quad(lambda t: α(t) * np.exp(-t*z), 0, 1000)
|
|||
|
|
#+end_src
|
|||
|
|
|
|||
|
|
#+RESULTS:
|
|||
|
|
|
|||
|
|
Now let's try calculating \(K_{11}\).
|
|||
|
|
First we define its laplace transform
|
|||
|
|
#+begin_src jupyter-python
|
|||
|
|
def K_l(z):
|
|||
|
|
return np.array([[z, Ω * np.ones_like(z)], [-(Ω + np.imag(α_l(z))), z]]) / (z ** 2 + Ω ** 2 + Ω * np.imag(α_l(z)))
|
|||
|
|
#+end_src
|
|||
|
|
|
|||
|
|
#+RESULTS:
|
|||
|
|
|
|||
|
|
Plotting it may help to give us some intuition.
|
|||
|
|
#+begin_src jupyter-python
|
|||
|
|
res = np.linspace(.1, 10, 100)
|
|||
|
|
ims = np.linspace(-10, 10, 100)
|
|||
|
|
zs = np.meshgrid(res, 1j* ims, indexing="ij")
|
|||
|
|
zs = (zs[0] + zs[1])
|
|||
|
|
image = plt.imshow(K_l(zs)[0,1].imag)
|
|||
|
|
plt.colorbar(image, ax=plt.gca(), location='right', anchor=(0, 0.3), shrink=0.7)
|
|||
|
|
#+end_src
|
|||
|
|
|
|||
|
|
#+RESULTS:
|
|||
|
|
:RESULTS:
|
|||
|
|
: <matplotlib.colorbar.Colorbar at 0x7fa8d2f60c70>
|
|||
|
|
[[file:./.ob-jupyter/ed423eb6eb8ae253acdfbf34ebc5c0318e43b0ca.svg]]
|
|||
|
|
:END:
|
|||
|
|
|
|||
|
|
|
|||
|
|
Let's do the inverse transform.
|
|||
|
|
#+begin_src jupyter-python
|
|||
|
|
import quadpy
|
|||
|
|
|
|||
|
|
def integrand(s, t):
|
|||
|
|
return 1 / (2j * np.pi) * K_l(s) * np.exp(s * t)
|
|||
|
|
|
|||
|
|
def γ(t):
|
|||
|
|
return 1j * t + 10 + 1.1 * t**2/np.sqrt(t**2 + 1)
|
|||
|
|
|
|||
|
|
def γ_dot(t):
|
|||
|
|
return 1.1 * (2*t/np.sqrt(t**2 + 1) - (t/np.sqrt(t**2 + 1))**3) + 1j
|
|||
|
|
|
|||
|
|
def K(t, lim):
|
|||
|
|
return quadpy.quad(
|
|||
|
|
lambda s: integrand(γ(s), t) * γ_dot(s),
|
|||
|
|
-lim,
|
|||
|
|
lim,
|
|||
|
|
epsabs=1e-3,
|
|||
|
|
epsrel=1e-3,
|
|||
|
|
)[0]
|
|||
|
|
|
|||
|
|
|
|||
|
|
t = np.linspace(0, 1, 100)
|
|||
|
|
#plt.plot(t, [K11(τ, 200, .1).real for τ in t])
|
|||
|
|
#plt.plot(t, [K11(τ, 200, 10).real for τ in t])
|
|||
|
|
#plt.plot(t, [K11(τ, 200, 20).real for τ in t])
|
|||
|
|
#plt.plot(t, [K11(τ, 400).real for τ in t])
|
|||
|
|
from scipy import interpolate
|
|||
|
|
K_i = interpolate.interp1d(t, [K(τ, 300) for τ in t], axis=0)
|
|||
|
|
#+end_src
|
|||
|
|
|
|||
|
|
#+RESULTS:
|
|||
|
|
|
|||
|
|
|
|||
|
|
#+begin_src jupyter-python
|
|||
|
|
t = np.linspace(0, 1, 1000)
|
|||
|
|
Ks = K_i(t).real
|
|||
|
|
for i, j in np.meshgrid([0, 1], [0, 1]):
|
|||
|
|
plt.plot(t, Ks[:, i, j])
|
|||
|
|
#+end_src
|
|||
|
|
|
|||
|
|
#+RESULTS:
|
|||
|
|
[[file:./.ob-jupyter/f30e618e0e7ac37663f69ca4fbebc83317e671ed.svg]]
|