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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
function dynamics!(
state,
prev_state,
params,
states,
plot_data,
)
# Unpack the parameters
N, K, beta, A, dt, num_steps = params
for j in 1:num_steps
next_state = zeros(N)
# update the left particle (set to rest)
state[1] = 0
# update the right particle (set to rest)
state[N] = 0
# update the middle particles
for i in 2:N-1
a = 2 * state[i] - prev_state[i]
b = dt * dt * K * (state[i+1] - 2 * state[i] + state[i-1])
c = dt * dt * beta * ((state[i+1] - state[i])^3 + (state[i-1] - state[i])^3)
this_mass = 1
if i == Int(N / 2)
# time = i * dt
# push!(plot_data, (state[i]))
this_mass = 1
end
next_state[i] = (a + b + c) / this_mass
end
# push!(states, next_state)
# update the previous state
for i in 1:N
prev_state[i] = state[i]
end
# update the current state
for i in 1:N
state[i] = next_state[i]
end
# if the middle mass is close to 2, print the time and position
if state[Int(N / 2)] > 1.995
println("Time: ", j * dt, " Positions: ", state[Int(N / 2)])
println("Other Positions: ", state, "\n")
push!(states, (j * dt, state))
end
if (j * dt % 100000) == 0
println("TIME UPDATE: ", j * dt)
end
end
end
function get_initial_state(
N,
modes,
beta,
A,
)
state = zeros(N)
# amp = A * sqrt(2 / (N - 1))
# for i in 2:N-1
# state[i] = amp * sin((modes * π * (i - 1)) / (N - 1))
# end
# set the middle to be the largest
state[Int(N / 2)] = 2
return state
end
function run_simulation(
N,
modes,
beta,
A,
dt,
num_steps,
plot_data,
)
states = []
state = get_initial_state(N, 1, beta, A)
prev_state = zeros(N)
for i in 1:N
prev_state[i] = state[i]
end
params = (N, modes, beta, A, dt, num_steps)
dynamics!(state, prev_state, params, states, plot_data)
return states
end
# Run the simulation
N = 10 # number of masses
beta = 0 # cubic string spring
K = 100 # spring constant
A = 10 # amplitude
modes = 3 # number of modes to plot
final_time = 10000000 # seconds
dt = 0.001 # seconds
num_steps = Int(final_time / dt)
params = (N, K, beta, A, dt, num_steps)
plot_data = []
s = run_simulation(N, K, beta, A, dt, num_steps, plot_data)
println("\nStates of interest:", s)
# plot the inital positions and the final positions
# using Plots
# plot(get_initial_state(N, 1, beta, A), label = "Initial", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", title = "First Three Modes")
# plot!(s[end], label = "Final", marker = :circle)
# plot the middle of the lattice over time
# time_steps = [i * dt for i in 1:num_steps]
# plot(time_steps, plot_data, label = "Middle Mass Over Time", xlabel = "Time", ylabel = "Displacement")
# savefig("t/plot_data.png")
# print the points close to the origional displacement
# output = []
# for i in 1:length(plot_data)
# if plot_data[i] > 1.975
# push!(output, i)
# println("Time: ", i * dt, " Position: ", plot_data[i])
# end
# end
# animate the s array of positions
# anim = @animate for i in 1:length(s)
# plot(s[i], label = "t = $(round(i * dt, digits = 3))", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", ylim = (-3, 3))
# end
# mp4(anim, "t/plz.mp4", fps = 30)
# println("Output: ", output)
|
