master-thesis/python/graveyard/richard_hops/test_derivative.org

• Setup
  • Jupyter
  • Matplotlib
  • Richard (old) HOPS
• Model Setup
  • BCF
    • Plot
  • Make stoch-procs

Setup

Jupyter

  %load_ext autoreload
  %autoreload 2
  %load_ext jupyter_spaces

Matplotlib

  import matplotlib
  import matplotlib.pyplot as plt

  #matplotlib.use("TkCairo", force=True)
  %gui tk
  %matplotlib inline
  plt.style.use('ggplot')

Richard (old) HOPS

  import hierarchyLib
  import hierarchyData
  import numpy as np

  from stocproc.stocproc import StocProc_FFT, StocProc_KLE
  import bcf
  from dataclasses import dataclass
  import scipy
  import scipy.misc
  import scipy.signal

Model Setup

Basic parameters.

  ω_c = 0  # center of spect. dens
  σ = 1  # breadth BCF
  t_max = 1
  seed = 100

BCF

   @dataclass
   class GaussBCF:
       σ: float
       ω_c: float

       def I(self, ω):
           return (
               np.exp(-(((ω - self.ω_c) / self.σ) ** 2) / 2)
               ,* 1
               / (np.sqrt(2 * np.pi) * self.σ)
           )

       def __call__(self, τ):
           return (
               np.exp(-(((τ - self.ω_c) / self.σ) ** 2) / 2)
               ,* 1
               / (np.sqrt(2 * np.pi) * self.σ)
           )

       def __bfkey__(self):
           return self.σ, self.ω_c


   @dataclass
   class GaussBCFD:
       α: GaussBCF

       def __call__(self, τ):
           # minus is important
           return self.α(τ)/self.α.σ**4 * (self.α.σ**2 - τ**2 + 2*τ*self.α.ω_c-self.α.ω_c**2)

       def __bfkey__(self):
           return self.α.__bfkey__()


  α = GaussBCF(σ, ω_c)
  α_td = GaussBCFD(α)

Plot

  %%space plot
  t = np.linspace(-t_max, t_max, 1000)
  ω = np.linspace(ω_c - 10, ω_c + 10, 1000)
  fig, axs = plt.subplots(2)
  axs[0].plot(t, np.real(α(t)))
  axs[0].plot(t, np.real(α_td(t)))
  #axs[1].plot(ω, α_td.I(ω))

<matplotlib.lines.Line2D	at	0x7f6440edd8e0>
<matplotlib.lines.Line2D	at	0x7f6440eddca0>

Make stoch-procs

Normal.

  η = StocProc_KLE(
      α,
      t_max,
      seed=seed,
      tol=1e-3
  )

  η_d = StocProc_KLE(
      α_td,
      t_max,
      seed=seed,
      tol=1e-3
  )

stocproc.method_kle - INFO - check 33 grid points
stocproc.method_kle - INFO - calc_ac 3.504%, fredholm 6.936%, integr_intp 1.506%, spline 6.846%, calc_diff 51.646%, rest 29.561%
stocproc.method_kle - INFO - auto ng SUCCESSFUL max diff 1.439e-04 < tol 1.000e-03 ng 33 num evec 3
stocproc.method_kle - INFO - check 33 grid points
stocproc.method_kle - INFO - calc_ac 9.515%, fredholm 4.745%, integr_intp 1.916%, spline 6.469%, calc_diff 49.757%, rest 27.598%
stocproc.method_kle - INFO - auto ng SUCCESSFUL max diff 7.083e-04 < tol 1.000e-03 ng 33 num evec 3
alpha_k is real
alpha_k is real

  η.new_process(seed=seed)
  η_d.new_process(seed=seed)
  #plt.plot(η.t, np.real(η()), label="orig")
  plt.plot(η.t, np.real(η_d()))
  plt.plot(η.t, scipy.misc.derivative(η, η.t, dx=1e-6).real, label="deriv")
  plt.legend()

stocproc.stocproc - INFO - use fixed seed (100) for new process
stocproc.stocproc - INFO - use fixed seed (100) for new process
<matplotlib.legend.Legend at 0x7f6440e64640>