-
Notifications
You must be signed in to change notification settings - Fork 2
/
rk4_demo.py
77 lines (58 loc) · 1.67 KB
/
rk4_demo.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
'''
Integration of a reaction-diffusion system exhibiting anomalous diffusion with a power-law memory kernel
This file initializes a set of parameters, runs the integration for a relatively generous amount of time, and then plots the output
William Gilpin, Spakowitz Group at Stanford University, 2015
'''
from matplotlib.pyplot import *
from scipy import *
from numpy import *
## pick parameter values
params = dict()
# width of reactive well
params['WELL_DIAM'] = .2
# width of overall potential well. smaller this is, the stronger the forcing
params['POT_DIAM'] = .05
params['KAPPA'] = 4e-1
params['DCOEFF'] = 1e-4
# alpha less than one
params['ALPHA']= .5
## set integrator settings
settings = dict()
space_pts = 25
ACTUAL_LENGTH = .5
dx = ACTUAL_LENGTH/space_pts
space = linspace(0.0, ACTUAL_LENGTH, space_pts)
space = space+dx
time_pts = 1e5
start_time = 0.0
stop_time = 17.0
dt = (stop_time-start_time)/time_pts
times = linspace(start_time, stop_time, time_pts)
times = times + dt
# initial conditions
y0 = ones(space_pts)
# settings['dx'] = dx
settings['space'] = space
settings['times'] = times
sol = rk4_mesh(y0, nxt_step, settings, params)
# Plot the results
# plot space-time concentration
figure()
imshow(sol, aspect='auto')
# plot time slices
figure()
solrange = len(sol[0,:])
num_slices = 50
hold(True)
slice_range = floor(solrange/num_slices)
plot(space, sol[:,::slice_range])
plot(space, sol[:,-1])
xlabel("Position along well")
ylabel("Probability")
# 'alpha_' + str(params['ALPHA']) + '__diff_' + str(params['DCOEFF'])
# savefig('harmonic_alpha_1__a_p2__V_p4.pdf')
# figure()
total_conc = sum(sol, axis = 0)
plot(times,total_conc[:-1])
ylabel("Total count")
xlabel("time")