notebook.community

Edit and run 





In [1]:



    


import numpy as np
import sympy as sp

%matplotlib nbagg
import matplotlib.pyplot as plt

from pyNT import Wiener, SDE, ODE



    











In [2]:



    


C1 = sp.Symbol('C_1')
C2 = sp.Symbol('C_2')
A  = sp.Symbol('A')



    











In [3]:



    


lambda1  = .2
lambda2  = .12
lambda21 = -.05
delta1   = .25
delta2   = .12
delta0   = .3

coeff_system = ODE(x = [C1, C2, A],
                   f = [-lambda1*C1 - lambda21*C2 + delta1,
                        -lambda2*C2 + delta2,
                        -(C1**2)/2 - (C2**2)/2 + delta0])



    











In [4]:



    


fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(111)
for solver in [['Heun', 'cRK'],
               ['cRK', 'Taylor4']]:
    evdt = coeff_system.get_evdt(
        X0 = np.zeros(3),
        solver = solver,
        h0 = .5,
        exp_range = range(10))
    # no errorbars since there's only one trajectory
    ax.plot(evdt[:, 0],
            evdt[:, 2],
            label = '{0} vs {1}'.format(solver[0], solver[1]))
ax.set_xscale('log')
ax.set_yscale('log')
ax.legend(loc = 'best')



    


















    




















    
















    
Out[4]:








<matplotlib.legend.Legend at 0x7f5cc1918630>

















In [5]:



    


coeff_time = np.linspace(.0, 1., 1025)
coeff_traj = coeff_system.cRK(2.**(-10), 2**10, np.zeros(3))



    











In [6]:



    


fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(111)
ax.plot(coeff_time, coeff_traj[:, 0], label = '$C_1$')
ax.plot(coeff_time, coeff_traj[:, 1], label = '$C_2$')
ax.plot(coeff_time, coeff_traj[:, 2], label = '$A$')
ax.legend(loc = 'best')



    


















    




















    
















    
Out[6]:








<matplotlib.legend.Legend at 0x7f5cc0b780b8>

















In [7]:



    


# reverse coefficients
traj_C1 = coeff_traj[::-1, 0]
traj_C2 = coeff_traj[::-1, 1]
traj_A  = coeff_traj[::-1, 2]



    











In [8]:



    


x = sp.Symbol('x')
y = sp.Symbol('y')
xy_sys = SDE(
    x = [x, y],
    a = [-lambda1*x,
         -lambda21*x - lambda2*y],
    b = [[1, 0],
         [0, 1]])



    











In [9]:



    


dt0, evdt0, deltae0 = xy_sys.get_evdt_vs_M(
    ntraj = 32,
    X0 = np.ones(2),
    h0 = 1.,
    #solver = ['EM', 'explicit_1p5_additive'],
    solver = ['EM', 'explicit_1p5_additive'],
    exp_range = range(7))



    


















    




















    
















    






<class 'int'>
<class 'int'>
<class 'int'>
<class 'int'>
<class 'int'>
<class 'int'>
<class 'int'>

















In [10]:



    


print(dt0)
print(evdt0)
print(deltae0)
print(2.**(-8))



    


















    






[0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25]
[0.00036563045439868267, 0.00072541226391117486, 0.0014317822856912126, 0.0029009909606799488, 0.0058797807662175064, 0.011813855032817747, 0.023912086668089118]
[1.637895962020072e-05, 3.3517110633991773e-05, 6.0168998747841968e-05, 0.0001344722960226638, 0.00026194844514277135, 0.00055297254297723283, 0.0010928776008797812]
0.00390625

















In [11]:



    


# get solution
nbatches = 20
ntraj = 64
# with dt = 2^(-8) = 0.0039..., the error is 0.0034 after integrating 1.
w = Wiener(nsteps = 2**8,
           dt = 2.**(-8),
           noise_dimension = xy_sys.get_noise_dimension(),
           solution_shape = [nbatches, ntraj])
w.initialize()
xy_traj = xy_sys.explicit_1p5_additive(w, np.ones(2))



    











In [12]:



    


# check that components of the Wiener process are different
print(w.W.shape)
print(w.dW.shape)
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(111)
ax.plot(w.get_time(), w.W[:, 0, 0, 0], label = '$X$')
ax.plot(w.get_time(), w.W[:, 1, 0, 0], label = '$Y$')
ax.legend(loc = 'best')



    


















    






(257, 2, 20, 64)
(256, 2, 20, 64)










    




















    
















    
Out[12]:








<matplotlib.legend.Legend at 0x7f5cc0a2cb70>

















In [13]:



    


print(xy_traj.shape)
traj_C1 = traj_C1[::4]
traj_C2 = traj_C2[::4]
traj_A  =  traj_A[::4]
print(traj_C1.shape)



    


















    






(257, 2, 20, 64)
(257,)

















In [14]:



    


X = xy_traj[:, 0]
Y = xy_traj[:, 1]



    











In [15]:



    


fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(111)
ax.plot(w.get_time(), X[:, 0, 0], label = '$X$')
ax.plot(w.get_time(), Y[:, 0, 0], label = '$Y$')
ax.legend(loc = 'best')
# the curves look very similar to the Wiener components,
#but they are different



    


















    




















    
















    
Out[15]:








<matplotlib.legend.Legend at 0x7f5cc09fcda0>

















In [16]:



    


# naive B
batch_pick = np.random.randint(0, nbatches)
traj_pick  = np.random.randint(0, ntraj)
B = np.exp(-X*traj_C1[:, None, None]
           -Y*traj_C2[:, None, None]
           -traj_A[:, None, None])
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(111)
ax.plot(w.get_time(), B[:, batch_pick, traj_pick], label = '$B$')
ax.legend(loc = 'best')



    


















    




















    
















    
Out[16]:








<matplotlib.legend.Legend at 0x7f5cc095e828>

















In [17]:



    


# compute error
# we're going to assume it explodes from 0
final_error = evdt0[0]
rel_error_plus = 1 + np.linspace(0,  final_error, X.shape[0])
rel_error_minus = 1 + np.linspace(0, -final_error, X.shape[0])
# note: I ignored the errors in C_1, C_2 and A since they're at least ten
#orders of magnitude less than the SDE errors...

fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(111)
ax.plot(w.get_time(), rel_error_plus)
ax.plot(w.get_time(), rel_error_minus)



    


















    




















    
















    
Out[17]:








[<matplotlib.lines.Line2D at 0x7f5cc09290b8>]

















In [18]:



    


batch_pick = np.random.randint(0, nbatches)
traj_pick  = np.random.randint(0, ntraj)
B = np.exp(-X*traj_C1[:, None, None]
           -Y*traj_C2[:, None, None]
           -traj_A[:, None, None])
Bpp = np.exp(
    -X*traj_C1[:, None, None]*rel_error_plus[:, None, None]
    -Y*traj_C2[:, None, None]*rel_error_plus[:, None, None]
    -traj_A[:, None, None])
Bpm = np.exp(
    -X*traj_C1[:, None, None]*rel_error_plus[:, None, None]
    -Y*traj_C2[:, None, None]*rel_error_minus[:, None, None]
    -traj_A[:, None, None])
Bmp = np.exp(
    -X*traj_C1[:, None, None]*rel_error_minus[:, None, None]
    -Y*traj_C2[:, None, None]*rel_error_plus[:, None, None]
    -traj_A[:, None, None])
Bmm = np.exp(
    -X*traj_C1[:, None, None]*rel_error_minus[:, None, None]
    -Y*traj_C2[:, None, None]*rel_error_minus[:, None, None]
    -traj_A[:, None, None])
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(111)
ax.plot(w.get_time(), B[:, batch_pick, traj_pick], label = '$B$')
ax.plot(w.get_time(), Bpp[:, batch_pick, traj_pick], label = '$B_{++}$')
ax.plot(w.get_time(), Bpm[:, batch_pick, traj_pick], label = '$B_{+-}$')
ax.plot(w.get_time(), Bmp[:, batch_pick, traj_pick], label = '$B_{-+}$')
ax.plot(w.get_time(), Bmm[:, batch_pick, traj_pick], label = '$B_{--}$')
ax.legend(loc = 'best')



    


















    




















    
















    
Out[18]:








<matplotlib.legend.Legend at 0x7f5cc0a42cc0>

Comment:

The result at $t=1$ is always 1 because we impose that $C_1$, $C_2$ and $A$ are 0 there (so we impose that $B(1) = 1$, it's independent of what methods we're using). The difference between different estimates of $B$ must therefore be 0 at $t=1$. It is also 0 at $t=0$ since the errors are 0 there as well.





In [19]:



    


# statistics on fluctuations in B
fullB = np.array([Bpp, Bpm, Bmp, Bmm])
minB = np.min(fullB, axis = 0)
maxB = np.max(fullB, axis = 0)
deltaB = maxB - minB
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(111)
ax.plot(w.get_time(), deltaB[:, batch_pick, traj_pick], label = '$\delta B$')
ax.legend(loc = 'best')



    


















    




















    
















    
Out[19]:








<matplotlib.legend.Legend at 0x7f5cc08d9e80>

















In [20]:



    


fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(111)
ax.plot(w.get_time(), np.average(deltaB, axis = (1, 2)), label = '$\delta B$')
ax.legend(loc = 'best')



    


















    




















    
















    
Out[20]:








<matplotlib.legend.Legend at 0x7f5cc0821470>

















In [21]:



    


xy_traj_alt = xy_sys.explicit_1p0(w, np.ones(2))



    











In [22]:



    


from pyNT.sde import get_t1ma_nm1
deltaB = np.abs(
    np.exp(
        -xy_traj[:, 0]*traj_C1[:, None, None]
        -xy_traj[:, 1]*traj_C2[:, None, None]
        -traj_A[:, None, None]) -
    np.exp(
        -xy_traj_alt[:, 0]*traj_C1[:, None, None]
        -xy_traj_alt[:, 1]*traj_C2[:, None, None]
        -traj_A[:, None, None]))
avdB = np.average(deltaB,
                  axis = (1, 2))
sigma2 = (np.sum((np.average(deltaB, axis = 2) - avdB[:, None])**2, axis = 1) /
          (deltaB.shape[1] - 1))
errB = (get_t1ma_nm1(0.99, deltaB.shape[1] - 1) *
        (sigma2 / deltaB.shape[1])**.5)
ax = plt.figure(figsize=(6,4)).add_subplot(111)
ax.plot(w.get_time(),
        avdB,
        label = '$\\langle \\delta B \\rangle$')
ax.plot(w.get_time(),
        np.max([np.min(deltaB,
                       axis = (1, 2)),
                avdB - errB], axis = 0),
        label = '$\\min \\delta B$')
ax.plot(w.get_time(),
        np.min([np.max(deltaB,
                       axis = (1, 2)),
                avdB + errB], axis = 0),
        label = '$\\max \\delta B$')
ax.legend(loc = 'best')



    


















    




















    
















    
Out[22]:








<matplotlib.legend.Legend at 0x7f5cc080b048>

Content source: chichilalescu/pyNT

Similar notebooks:

notebook.community | gallery | about