hiro/stocproc
1
0
Fork 0
mirror of https://github.com/vale981/stocproc synced 2025-03-04 17:21:42 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

made StocProc_TS more robust for w=0 singular SD

Browse source
This commit is contained in:
Richard Hartmann 2019-01-19 12:20:33 +01:00
parent 4054867f34
commit d078ac32e1
2 changed files with 38 additions and 8 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

32
stocproc/method_ft.py
View file

@ -140,17 +140,41 @@ def yk(h, k):
# I, feed_back = fourier_integral_TanhSinh_with_feedback(integrand, w_max, tau, h, kmax)
# return I
def get_x_w_and_dt(n, x_max):
t_l, d_t = np.linspace(-3, 3, n, retstep=True)
def get_t_max_for_singularity_ts(f, a, b, tol):
"""
chose tmax such that |w_(tmax) I(g(tmax))| < tol
"""
sc = (b-a)/2
t_max = 3
while t_max < 6:
s_tmax = np.sinh(t_max) * np.pi / 2
g_tmax = 1 / (np.exp(s_tmax) * np.cosh(s_tmax))
w_tmax = np.pi * np.cosh(t_max) / 2 / np.cosh(s_tmax) ** 2
f_x = f(a + sc*g_tmax)
tmp = abs(sc * f_x * w_tmax)
if tmp < tol:
#print("for t_max {} (boundary at singulatigy) error condition fulfilled (err est {} < tol {})".format(t_max, tmp, tol))
return t_max
else:
#print("for t_max {} (boundary at singulatigy) got err est {} >= tol {} -> increase t_max".format(t_max, tmp, tol))
pass
t_max += 0.5
def get_x_w_and_dt(n, x_max, t_max):
t_l, d_t = np.linspace(-3, t_max, n, retstep=True)
s_t = np.sinh(t_l) * np.pi / 2
x = x_max / 2 / (np.exp(s_t) * np.cosh(s_t))
w = x_max / 2 * d_t * np.pi * np.cosh(t_l) / 2 / np.cosh(s_t) ** 2
return x, w
def fourier_integral_TanhSinh(f, x_max, n, tau_l):
def fourier_integral_TanhSinh(f, x_max, n, tau_l, t_max_ts):
_intgr = lambda x, tau: f(x)*np.exp(-1j*x*tau)
x, w = get_x_w_and_dt(n, x_max)
x, w = get_x_w_and_dt(n, x_max, t_max_ts)
I = np.asarray([np.sum(_intgr(x, tau) * w) for tau in tau_l])
return I

14
stocproc/stocproc.py
View file

@ -486,9 +486,12 @@ class StocProc_TanhSinh(_abcStocProc):
log.info("yields N = {} (time domain)".format(N))
log.info("find accurate discretisation in frequency domain")
wmax = method_ft.find_integral_boundary(spectral_density, tol=intgr_tol/4, ref_val=1, max_val=1e6, x0=0.777)
wmax = method_ft.find_integral_boundary(spectral_density, tol=intgr_tol/10, ref_val=1, max_val=1e6, x0=0.777)
log.info("wmax:{}".format(wmax))
t_max_ts = method_ft.get_t_max_for_singularity_ts(lambda w: spectral_density(w)/np.pi,
0, wmax, intgr_tol)
tau = np.linspace(0, t_max, 35)
n = 16
d = intgr_tol+1
@ -497,8 +500,10 @@ class StocProc_TanhSinh(_abcStocProc):
I = method_ft.fourier_integral_TanhSinh(f = lambda w: spectral_density(w)/np.pi,
x_max=wmax,
n=n,
tau_l=tau)
tau_l=tau,
t_max_ts = t_max_ts)
bcf_ref_t = alpha(tau)
d = np.abs(bcf_ref_t-I)/abs(bcf_ref_t[0])
d = np.max(d)
print("n:{} d:{} tol:{}".format(n, d, intgr_tol))
@ -508,7 +513,8 @@ class StocProc_TanhSinh(_abcStocProc):
num_FT = method_ft.fourier_integral_TanhSinh(f = lambda w: spectral_density(w)/np.pi,
x_max=wmax,
n=n,
tau_l=tau)
tau_l=tau,
t_max_ts=t_max_ts)
bcf_ref_t = alpha(tau)
d = np.max(np.abs(num_FT - bcf_ref_t) / np.abs(bcf_ref_t[0]))
@ -529,7 +535,7 @@ class StocProc_TanhSinh(_abcStocProc):
plt.show()
assert d <= intgr_tol, "d:{}, intgr_tol:{}".format(d, intgr_tol)
yk, wk = method_ft.get_x_w_and_dt(n, wmax)
yk, wk = method_ft.get_x_w_and_dt(n, wmax, t_max_ts)
self.omega_k = yk
self.fl = np.sqrt(wk*spectral_density(self.omega_k)/np.pi)