In [2]:

    

# Libraries are in parent directory
import sys
sys.path.append('../')

from dtrw import *
import numpy as np
import matplotlib.pyplot as plt


    

In [4]:

    

alpha_1 = 0.5
D_alpha = 1.
T = 0.1
k_1 = 1.

dXs = [0.1, 0.05, 0.02, 0.01]
dTs = [.001, .0001, .00001, .000001]

dTs = [1e-5, 1e-5, 1e-6, 1e-6]
nXs = [5, 10, 20, 40, 80]#[21, 41, 61, 81]
dXs = [0.5, 0.2, 1.0/7., 0.1, 1.0/15, 0.05, 0.025, 0.0125]
#for dT in dTs:

for j in range(len(dXs)):
    
    #nX  = nXs[j]
    #dT = dTs[j]
    #dX = math.sqrt(2. * D_alpha * pow(dT, alpha_1))
    
    #dX = 2. / nX
    dX = dXs[j]
    dT = pow(dX * dX / (2.0 * D_alpha), 1.0 / alpha_1)
    print "dT:", dT
    print "dX:", dX
    
    L = math.floor(1.0 / dX)
    print L, 1.0 / dX
    print round(T / dT), T / dT
    print ""
    
    r = 2. * D_alpha * pow(dT, alpha_1) / (dX * dX)


    

    




dT: 0.015625
dX: 0.5
2.0 2.0
6.0 6.4

dT: 0.0004
dX: 0.2
5.0 5.0
250.0 250.0

dT: 0.000104123281966
dX: 0.142857142857
7.0 7.0
960.0 960.4

dT: 2.5e-05
dX: 0.1
10.0 10.0
4000.0 4000.0

dT: 4.93827160494e-06
dX: 0.0666666666667
15.0 15.0
20250.0 20250.0

dT: 1.5625e-06
dX: 0.05
20.0 20.0
64000.0 64000.0

dT: 9.765625e-08
dX: 0.025
40.0 40.0
1024000.0 1024000.0

dT: 6.103515625e-09
dX: 0.0125
80.0 80.0
16384000.0 16384000.0

In [26]:

    

dT * 1024 * 1000


    

    
Out[26]:





0.10000000000000005

In [3]:

    

I_sol_21 = np.loadtxt("Solution_alpha0p8_k1_dx2ov21_T0p1.csv", delimiter=",")
I_sol_41 = np.loadtxt("Solution_alpha0p8_k1_dx2ov41_T0p1.csv", delimiter=",")
I_sol_61 = np.loadtxt("Solution_alpha0p8_k1_dx2ov61_T0p1.csv", delimiter=",")
I_sol_81 = np.loadtxt("Solution_alpha0p8_k1_dx2ov81_T0p1.csv", delimiter=",")


    

In [4]:

    

a_21 = np.loadtxt("a_DTRW_dT_0.001000_dX_0.095238_alpha_0.800000.csv", delimiter=",")
b_21 = np.loadtxt("b_DTRW_dT_0.001000_dX_0.095238_alpha_0.800000.csv", delimiter=",")
a_41 = np.loadtxt("a_DTRW_dT_0.000100_dX_0.048780_alpha_0.800000.csv", delimiter=",")
b_41 = np.loadtxt("b_DTRW_dT_0.000100_dX_0.048780_alpha_0.800000.csv", delimiter=",")
a_61 = np.loadtxt("a_DTRW_dT_0.000050_dX_0.032787_alpha_0.800000.csv", delimiter=",")
b_61 = np.loadtxt("b_DTRW_dT_0.000050_dX_0.032787_alpha_0.800000.csv", delimiter=",")
a_81 = np.loadtxt("a_DTRW_dT_0.000010_dX_0.024691_alpha_0.800000.csv", delimiter=",")
b_81 = np.loadtxt("b_DTRW_dT_0.000010_dX_0.024691_alpha_0.800000.csv", delimiter=",")


    

In [12]:

    

I_trev = np.loadtxt("trev_cable_eq_results/Soln_at_t01.dat")

a_10 = np.loadtxt("a_DTRW_dT_0.000010_dX_0.200000_alpha_0.500000.csv", delimiter=",")
b_10 = np.loadtxt("b_DTRW_dT_0.000010_dX_0.200000_alpha_0.500000.csv", delimiter=",")
a_20 = np.loadtxt("a_DTRW_dT_0.000010_dX_0.100000_alpha_0.500000.csv", delimiter=",")
b_20 = np.loadtxt("b_DTRW_dT_0.000010_dX_0.100000_alpha_0.500000.csv", delimiter=",")
a_40 = np.loadtxt("a_DTRW_dT_0.000001_dX_0.050000_alpha_0.500000.csv", delimiter=",")
b_40 = np.loadtxt("b_DTRW_dT_0.000001_dX_0.050000_alpha_0.500000.csv", delimiter=",")
    
print a_10.shape, np.linspace(-1.,1.,11, endpoint=True).shape
print a_20.shape, np.linspace(-1.,1.,21, endpoint=True).shape
print a_40.shape, np.linspace(-1.,1.,41, endpoint=True).shape

xs = np.arange(-1., 1.05, 0.05)
plt.plot(xs, I_trev[:,1], label='Analytic soln')
lin1 = plt.plot(np.linspace(-1.,1.,11, endpoint=True), (a_10+b_10), 'x', label='DTRW dX=1/5')
lin1 = plt.plot(np.linspace(-1.,1.,11, endpoint=True), b_10, 'x', label='DTRW dX=1/5')
lin2 = plt.plot(np.linspace(-1.,1.,21, endpoint=True), (a_20+b_20), 'v', label='DTRW dX=1/10')
lin3 = plt.plot(np.linspace(-1.,1.,41, endpoint=True), (a_40+b_40), 'o', label='DTRW dX=1/20')
plt.legend()
plt.xlim([-1.,1.])
plt.show()


    

    




(11,) (11,)
(21,) (21,)
(41,) (41,)

In [9]:

    

I_trev = np.loadtxt("trev_cable_eq_results/Soln_at_t01.dat")

a_10 = np.loadtxt("a_DTRW_dT_0.000010_dX_0.200000_alpha_0.500000.csv", delimiter=",")
b_10 = np.loadtxt("b_DTRW_dT_0.000010_dX_0.200000_alpha_0.500000.csv", delimiter=",")
a_20 = np.loadtxt("a_DTRW_dT_0.000010_dX_0.100000_alpha_0.500000.csv", delimiter=",")
b_20 = np.loadtxt("b_DTRW_dT_0.000010_dX_0.100000_alpha_0.500000.csv", delimiter=",")
a_40 = np.loadtxt("a_DTRW_dT_0.000001_dX_0.050000_alpha_0.500000.csv", delimiter=",")
b_40 = np.loadtxt("b_DTRW_dT_0.000001_dX_0.050000_alpha_0.500000.csv", delimiter=",")

print a_40.shape


    

    




(41,)

In [16]:

    

plt.plot(I_sol_21 * np.exp(-0.1))
plt.plot((a_21+b_21)[:,-1])
plt.plot(I_sol_41 * np.exp(-0.1))
plt.plot((a_41+b_41)[:,-1])
plt.plot(I_sol_61 * np.exp(-0.1))
plt.plot((a_61+b_61)[:,-1])
plt.plot(I_sol_81 * np.exp(-0.1))
plt.plot((a_81+b_81)[:,-1])
plt.show()


    

In [9]:

    

print I_trev[:,1] / (a_41+b_41)[:,-1]


    

    




[ 1.27761439  1.26483523  1.24363341  1.2155739   1.18259474  1.14672723
  1.10986556  1.07362026  1.03925699  1.00770148  0.97958455  0.95530487
  0.93509347  0.91907175  0.9072994   0.8998123   0.89665187  0.8978879
  0.90363721  0.91408011  0.92947659  0.91408011  0.90363721  0.8978879
  0.89665187  0.8998123   0.9072994   0.91907175  0.93509347  0.95530487
  0.97958455  1.00770148  1.03925699  1.07362026  1.10986556  1.14672723
  1.18259474  1.2155739   1.24363341  1.26483524  1.2776144 ]

In [ ]: