-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbastantes_osci.jl~
98 lines (82 loc) · 3.17 KB
/
bastantes_osci.jl~
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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
using OrdinaryDiffEq, ForwardDiff, Distributions, RandomNumbers, NPZ
function initial_sampling(Eb,w1,w2,mb)
nums=rand(1)*2*pi
x1=cos(nums[1])*sqrt(Eb/(0.5*mb*w1*w1))
y1=sin(nums[1])*sqrt(2*mb*Eb)
x2=cos(nums[1])*sqrt(Eb/(0.5*mb*w2*w2))
y2=sin(nums[1])*sqrt(2*mb*Eb)
return x1,y1,-x2,-y2
end
function bolz_energy(a)
return rand(Truncated(Exponential(1), 0., 10))
end
function bolz_freq(a)
return rand()*0.4324555320336759+0.63
end
initial_position = [0.1, 0.,0.,0.,0.,0.,0.,0.,0.]
initial_momentum=[0. ,0.1,0.,0.,0.,0.,0.,0.,0.]
#numeros=[6,8,10,16,20,26,30,36,40,50]
numeros=[36,40,50]
for ii in 1:size(numeros)[1]
println(ii)
ii=Int(ii)
N=numeros[ii]
initial_position=ones(N+1)
initial_momentum=ones(N+1)
a=0.2
b=0.01
m=[1.,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]
w=[0.,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8]
g=[0.,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]
X0=17.
P0=1.
golden=(1+sqrt(big(5)))/2
Es=zeros(N+1)
time=100000.
#time=1.
tspan=(0., time)
#fig=figure()
w=zeros(N+1)
g=zeros(N+1)
E=zeros(N+1)
m=ones(N+1) .*0.1
m[1]=1.
w=w .* 0. .+ 1.
w=(bolz_freq.(w))
Es=(bolz_energy.(Es))
for i in 1:(N-1)/2
i=Int(i)
w[i*2+1]=w[i*2]+0.0001
end
#println(w)
W=maximum(w)
g=g .* 0. .+ 0.3/N
dts=0.01*2. *pi/W
#Hsys(q,p) = p[1]^2. /(2. *m[1]) +m[1]*w[1]^2. *q[1]^2. /2. +0.5* -q[1]*g[2]*q[2]+p[2]^2. /(2. *m[2])+m[2]*w[2]^2. *q[2]^2. /2. + q[1]^2. *g[2]^2. /(2. *m[2]*w[2]^2.) -q[1]*g[3]*q[3]+p[3]^2. /(2. *m[3])+m[3]*w[3]^2. *q[3]^2. /2. + q[1]^2. *g[3]^2. /(2. *m[3]*w[3]^2.) -q[1]*g[4]*q[4]+p[4]^2. /(2. *m[4])+m[4]*w[4]^2. *q[4]^2. /2. + q[1]^2. *g[4]^2. /(2. *m[4]*w[4]^2.) -q[1]*g[5]*q[5]+p[5]^2. /(2. *m[5])+m[5]*w[5]^2. *q[5]^2. /2. + q[1]^2. *g[5]^2. /(2. *m[5]*w[5]^2.) -q[1]*g[6]*q[6]+p[6]^2. /(2. *m[6])+m[6]*w[6]^2. *q[6]^2. /2. + q[1]^2. *g[6]^2. /(2. *m[6]*w[6]^2.) -q[1]*g[7]*q[7]+p[7]^2. /(2. *m[7])+m[7]*w[7]^2. *q[7]^2. /2. + q[1]^2. *g[7]^2. /(2. *m[7]*w[7]^2.) -q[1]*g[8]*q[8]+p[8]^2. /(2. *m[8])+m[8]*w[8]^2. *q[8]^2. /2. + q[1]^2. *g[8]^2. /(2. *m[8]*w[8]^2.)-q[1]*g[9]*q[9]+p[9]^2. /(2. *m[9])+m[9]*w[9]^2. *q[9]^2. /2. + q[1]^2. *g[9]^2. /(2. *m[9]*w[9]^2.)
function Hsys(q,p,N0)
result=p[1]^2. /(2. *m[1]) - a*q[1]^2. /2. +b*q[1]^4. /4.
for i in 1:N0
result += -q[1]*g[i+1]*q[i+1]+p[i+1]^2. /(2. *m[2])+m[i+1]*w[i+1]^2. *q[i+1]^2. /2. + q[1]^2. *g[i+1]^2. /(2. *m[i+1]*w[i+1]^2.)
end
return result
end
pdot(dp,p,q,params,t) = ForwardDiff.gradient!(dp, q->-Hsys(q, p,N), q)
qdot(dq,p,q,params,t) = ForwardDiff.gradient!(dq, p-> Hsys(q, p,N), p)
initial_position=initial_position .* 0.
initial_momentum=initial_momentum .* 0.
initial_position[1]=0.0 #xz0[jj]
initial_momentum[1]=0.0 #pz1[jj]
for i in 1:(N-1)/2
i=Int(i)
initial_position[i*2], initial_momentum[i*2],initial_position[i*2+1], initial_momentum[i*2+1]=initial_sampling(Es[i*2],w[i*2],w[i*2+1],m[i*2])
end
prob = DynamicalODEProblem(pdot, qdot, initial_momentum, initial_position, tspan)
@time sol = solve(prob, CalvoSanz4(), dt=dts)
#subplot(2, 1, 1)
#plot(sol[N+2,:],sol[1,:],color="b",linewidth=0.5,alpha=0.8)
#subplot(2, 1, 2)
#plot(sol.t[:], sol[N+2,:], color="b",linewidth=0.5,alpha=0.8)
npzwrite("x$N.npy", sol[N+2,:])
npzwrite("p$N.npy", sol[1,:])
npzwrite("t$N.npy", sol.t[:])
end