diff options
Diffstat (limited to 't/plz.jl')
| -rw-r--r-- | t/plz.jl | 767 |
1 files changed, 0 insertions, 767 deletions
diff --git a/t/plz.jl b/t/plz.jl deleted file mode 100644 index 962398a..0000000 --- a/t/plz.jl +++ /dev/null @@ -1,767 +0,0 @@ -function calculate_force( - left_pos, - middle_pos, - right_pos, - K, - alpha = 0.0, - beta = 0.0, -) - linear_force = K * (left_pos + right_pos - 2 * middle_pos) - quadratic_force = alpha * (left_pos - middle_pos)^2 + alpha * (right_pos - middle_pos)^2 - cubic_force = beta * (left_pos - middle_pos)^3 + beta * (right_pos - middle_pos)^3 - - return linear_force + quadratic_force + cubic_force -end - -function tendency!(du, u, params, t) - # unpack the params - N, modes, beta, A, dt, num_steps = params - - # get the positions and momenta - qs = u[1:2:end] - ps = u[2:2:end] - - num_masses = length(qs) - - # go over the points in the lattice and update the state - for i in 2:num_masses-1 - # left_index = max(1, i - 1) - # right_index = min(N, i + 1) - - du[i*2-1] = ps[i] - force = calculate_force(qs[i-1], qs[i], qs[i+1], 1, 0, beta) - du[i*2] = force - end - - # make last point same as first at rest - du[num_masses*2-1] = du[1] = 0 # set to 0 - du[num_masses*2] = du[2] = 0 - - # if ps[Int(N / 2)] / M >= 1 - # println("(in sim!) Time: ", t, " Vel: ", ps[Int(N / 2)] / M) - # # println("Other Positions: ", qs) - # println("Other Velocities: ", ps, "\n") - # end -end - -function get_initial_state( - N, - modes, - beta, - A, -) - state = zeros(N + 2) - amp = A * sqrt(2 / (N + 1)) - for i in 2:N+1 - state[i] = amp * sin((modes * π * (i - 1)) / (N + 1)) - end - return state -end - -using DifferentialEquations -function run_simulation( - N, - modes, - beta, - A, - dt, - num_steps, -) - println("Running simulation with params: N=$N, modes=$modes, beta=$beta, A=$A, dt=$dt, num_steps=$num_steps") - # 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) - - # state above is position, need to add in momentum - s_0 = zeros(2 * (N + 2)) - for i in 1:length(state) - s_0[i*2-1] = state[i] - end - final_time = num_steps * dt - t_span = (0.0, final_time) - prob = ODEProblem(tendency!, s_0, t_span, params) - sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8, maxiters = 1e10) # control simulation - - return sol -end - -function get_pos_from_state(state) - return state[1:2:end] -end - -# Run the simulation -# N = 32 # number of masses -# beta = 2.3 # cubic string spring -# A = 10 # amplitude -# modes = 3 # number of modes to plot -# final_time = 100000 # seconds -# dt = 0.05 # seconds -# num_steps = Int(final_time / dt) -# params = (N, modes, beta, A, dt, num_steps) -# println("Running simulation with params: ", params) -# sol = run_simulation(N, modes, beta, A, dt, num_steps) -# states = sol.u -# timesteps = sol.t -# println("Final state: ", states[end]) - -# plot the inital positions and the final positions -using Plots -# init_pos = get_initial_state(N, 1, beta, A) -# final_pos = get_pos_from_state(states[end]) -# println("Initial: ", init_pos) -# println("Final: ", final_pos) -# plot(init_pos, label = "Initial", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", title = "First Three Modes") -# plot!(final_pos, label = "Final", marker = :circle) -# println("Saving initial-final-beta15.png") -# savefig("t/initial-final-beta-$beta.png") - -function animate_states(states, timesteps, from = 1, to = 1000) - println("\nAnimating states") - # animate the s array of positions - # find the timestep that maps to the state index - i_start = 1 - i_end = length(timesteps) - for i in 1:length(timesteps) - if timesteps[i] <= from - i_start = i - continue - end - if timesteps[i] >= to - i_end = i - break - end - end - - anim = @animate for i in i_start:i_end - s = get_pos_from_state(states[i]) - t = timesteps[i] - plot(s, label = "t = $t", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", ylim = (-3, 3)) - end - - mp4(anim, "t/animated-fpu.mp4", fps = 30) - println("Done animating, wrote to animated-fpu.mp4\n") -end -# animate_states(states, timesteps, 149500, 150000) - -function caluclate_energies_for_mode(states, mode) - N = length(states[1]) / 2 - 2 - - total = [] - kinetic = [] - potential = [] - - for i in 1:length(states) - total_energy = 0 - kinetic_energy = 0 - potential_energy = 0 - - # find a_k for this mode - poses = get_pos_from_state(states[i]) - a_k = 0 - for j in 1:length(poses) - a_k += poses[j] * sin((mode * (j - 1) * π) / (N + 1)) - end - amp = sqrt(2 / (N + 1)) - a_k *= amp - - # get the potenital energy - omega_k = 2 * sin((mode * π) / (2 * (N + 1))) - p = 0.5 * omega_k^2 * a_k^2 - - # find the kinetic energy - v_k = 0 - vels = states[i][2:2:end] - for j in 1:length(vels) - v_k += vels[j] * sin((mode * (j - 1) * π) / (N + 1)) - end - v_k *= amp - k = 0.5 * v_k^2 - - - total_energy += k + p - kinetic_energy += k - potential_energy += p - push!(total, total_energy) - push!(kinetic, kinetic_energy) - push!(potential, potential_energy) - end - return total, kinetic, potential -end - -function calculate_total_energy_by_state(states, timesteps) - total = [] - kinetic = [] - potential = [] - for i in 1:length(states) - vels = states[i][2:2:end] - k = 0.5 * sum(vels .^ 2) - - poses = get_pos_from_state(states[i]) - p = 0 - for j in 2:length(poses) - p += beta * (poses[j] - poses[j-1])^2 - end - - push!(total, k + p) - push!(kinetic, k) - push!(potential, p) - end - - # plot all three over time on the same plot - plot(timesteps, total, label = "Total Energy", xlabel = "Time", ylabel = "Energy", title = "Total Energy Over Time") - plot!(timesteps, kinetic, label = "Kinetic Energy") - plot!(timesteps, potential, label = "Potential Energy") - - return total, kinetic, potential -end - -function print_energies_for_modes(states, max_mode = 3, index = 1) - # get the sum of total energy over all mode - total = 0 - kinetic = 0 - potential = 0 - for i in 1:max_mode - t, k, p = caluclate_energies_for_mode(states, i) - total += t[index] - kinetic += k[index] - potential += p[index] - end - println("Total Energy: ", total) - println("Kinetic Energy: ", kinetic) - println("Potential Energy: ", potential) -end - -function plot_phase_space(states, mass_number, intersecting_mass = 4) - # get the positions and momenta - qs = [] - ps = [] - - for i in 2:length(states) - # only store if the other mass is crossing 0 - if (states[i-1][intersecting_mass*2-1] < 0 && states[i][intersecting_mass*2-1] > 0) || (states[i-1][intersecting_mass*2-1] > 0 && states[i][intersecting_mass*2-1] < 0) - s = states[i] - q = s[mass_number*2-1] - p = s[mass_number*2] - push!(qs, q) - push!(ps, p) - end - end - - s = scatter( - qs, - ps, - label = "Mass $mass_number", - xlabel = "Displacement", - ylabel = "Momentum", - title = "Phase Space for Masses when Mass $intersecting_mass crosses 0", - legend = :topleft, - marker = :circle, - markersize = 2, - xlim = (-1.5, 1.5), - ylim = (-1.5, 1.5), - ) - - # savefig(s, "t/phase-space-mass$mass_number-beta15.png") - - return s -end - -function analyze_energies_of_n_modes(modes::Array{Int}, states, timesteps, beta) - println("Analyzing energies of modes: ", modes) - - total_by_mode = Dict() - kinetic_by_mode = Dict() - potential_by_mode = Dict() - for j in modes - t, k, p = caluclate_energies_for_mode(states, j) - - total_by_mode[j] = t - kinetic_by_mode[j] = k - potential_by_mode[j] = p - end - - total_energy = [] - kinetic_energy = [] - potential_energy = [] - for i in 1:length(states) - total = 0 - kinetic = 0 - potential = 0 - - for j in modes - total += total_by_mode[j][i] - kinetic += kinetic_by_mode[j][i] - potential += potential_by_mode[j][i] - end - - push!(total_energy, total) - push!(kinetic_energy, kinetic) - push!(potential_energy, potential) - end - - # plot all three over time on the same plot - plot(timesteps, total_energy, label = "Total Energy", xlabel = "Time", ylabel = "Energy", title = "Total Energy Over Time") - plot!(timesteps, kinetic_energy, label = "Kinetic Energy") - plot!(timesteps, potential_energy, label = "Potential Energy") - savefig("t/procs/energies-modes-sum-beta-$beta.png") - println("Saved energies-modes-sum-beta-$beta.png. Summary Below:\n") - - # plot the energies for each mode over time - plt = plot(title = "Energy Over Time for Modes $(join(modes, ", "))", xlabel = "Time", ylabel = "Energy") - for mode in modes - plot!(plt, timesteps, total_by_mode[mode], label = "Mode $mode") - # plot!(plt, timesteps, kinetic_by_mode[mode], label = "Mode $mode Kinetic") - # plot!(plt, timesteps, potential_by_mode[mode], label = "Mode $mode Potential") - end - savefig(plt, "t/procs/energies-modes-separate-beta-$beta.png") - - # print the energies - println("Total Energy Initial: ", total_energy[1]) - println("Total Energy Final: ", total_energy[end]) - println("Kinetic Energy Initial: ", kinetic_energy[1]) - println("Kinetic Energy Final: ", kinetic_energy[end]) - println("Potential Energy Initial: ", potential_energy[1]) - println("Potential Energy Final: ", potential_energy[end], "\n") - - println("Done analyzing energies of modes: ", modes, "\n") - - return total_energy, kinetic_energy, potential_energy, total_by_mode, kinetic_by_mode, potential_by_mode -end - -function animate_masses_phase_space(states, mass_numbers, intersecting_mass = 4) - println("\nAnimating phase space for mass $mass_numbers when mass $intersecting_mass crosses 0") - # animate the s array of positions - # find the timestep that maps to the state index - anim = @animate for i in mass_numbers - s = plot_phase_space(states, i, intersecting_mass) - s - end - - mp4(anim, "t/animated-phase-space-masses-beta-$beta.mp4", fps = 24) -end - -function plot_mode_phase_space(states, mode_num, intersecting_mode, beta) - println("\nPlotting phase space for mode $mode_num") - - # convert the state to a mode - a_modes = [] - vel_modes = [] - a_inter_modes = [] - vel_inter_modes = [] - - N = length(states[1]) / 2 - 2 - - for i in 1:length(states) - a_mode = 0 - vel_mode = 0 - a_inter_mode = 0 - vel_inter_mode = 0 - - positions = states[i][1:2:end] - velocities = states[i][2:2:end] - for j in 2:length(positions)-1 - a_mode += positions[j] * sin((mode_num * (j - 1) * π) / (N + 1)) - vel_mode += velocities[j] * sin((mode_num * (j - 1) * π) / (N + 1)) - a_inter_mode += positions[j] * sin((intersecting_mode * (j - 1) * π) / (N + 1)) - vel_inter_mode += velocities[j] * sin((intersecting_mode * (j - 1) * π) / (N + 1)) - end - - amp = sqrt(2 / (N + 1)) - a_mode *= amp - vel_mode *= amp - a_inter_mode *= amp - vel_inter_mode *= amp - - if i == 1 - push!(a_inter_modes, a_inter_mode) - push!(vel_inter_modes, vel_inter_mode) - continue - end - - if ((a_inter_modes[end] < 0 && a_inter_mode > 0) || (a_inter_modes[end] > 0 && a_inter_mode < 0)) - push!(a_modes, a_mode) - push!(vel_modes, vel_mode) - end - push!(a_inter_modes, a_inter_mode) - push!(vel_inter_modes, vel_inter_mode) - end - - - # plot the phase space - scatter(a_modes, vel_modes, label = "Beta = $beta", xlabel = "Displacement", ylabel = "Momentum", - title = "Phase Space for Mode $mode_num when a_$intersecting_mode crosses 0", - xlim = (-10, 10), - ylim = (-1, 1), - color = :red, - legend = :topright, - ) - # savefig("t/frames/phase-space-mode-$mode_num-beta-$beta.png") - println("Saved phase-space-mode-$mode_num-beta-$beta.png") -end - -function animate_phase_space_plots_over_beta(betas, t_max, desired_mode, intersecting_mode) - N = 32 # number of masses - A = 10 # amplitude - modes = 3 # number of modes to plot - final_time = t_max # seconds - dt = 0.05 # seconds - num_steps = Int(final_time / dt) - - results = [] - for b in betas - s = run_simulation(N, modes, b, A, dt, num_steps) - states = s.u - timesteps = s.t - - push!(results, states) - - println("Simualted beta: ", b, " to t=", timesteps[end]) - end - - anim = @animate for i in 1:length(betas) - plot_mode_phase_space(results[i], desired_mode, intersecting_mode, betas[i]) - end - - mp4(anim, "t/animated-phase-space-modes-over-beta.mp4", fps = 15) -end - -# animate_phase_space_plots_over_beta(collect(0:0.025:4), 25000, 1, 3) - -# analyze_energies_of_n_modes([1, 3, 5, 7, 9, 11], timesteps) - -# animate_masses_phase_space(states, collect(2:N+1), 17) - -# plot_mode_phase_space(states, 1, 3) - -function estimate_lynapov_exponent(t_max, beta, jitter = 0.001) - println("Estimating Lynapov Exponent for beta: ", beta) - - N = 32 # number of masses - A = 10 # amplitude - modes = 3 # number of modes to plot - println("t_msx: ", t_max) - final_time = t_max # seconds - dt = 0.05 # seconds - num_steps = Int(final_time / dt) - - t_span = (0.0, final_time) - params = (N, modes, beta, A, dt, num_steps) - - # run a simulation with normal parameters - is = get_initial_state(N, 1, beta, A) - # plot iniital state - plot(is, label = "Initial State", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", title = "Initial Positions of Masses", markersize = 4) - is_0 = zeros(2 * (N + 2)) - for i in 1:length(is) - is_0[i*2-1] = is[i] - end - prob_s = ODEProblem(tendency!, is_0, t_span, params) - # solve so that the timestep is exactly set to dt - s = solve(prob_s, SSPRK33(), dt = dt) - states_control = s.u - timesteps = s.t - - - # run a simulation with a small perturbation - # dynamics!(state, prev_state, params, states) - # state above is position, need to add in momentum - is2 = get_initial_state(N, 1, beta, A) .* (1 - jitter) # perturb the initial state - # is2 = get_initial_state(N, 1, beta, A) .- jitter # perturb the initial state - plot!(is2, label = "Perturbed Initial State", marker = :cross, xlabel = "Mass Number", ylabel = "Displacement", title = "Perturbed Initial Positions of Masses", msw = 2) - savefig("t/new/inits-jitter.png") - s_0 = zeros(2 * (N + 2)) - for i in 1:length(is2) - s_0[i*2-1] = is2[i] - end - - # add random jitter to the initial state - # middle_index = Int(N / 2) * 2 - # s_0[middle_index] -= jitter # subtract a little from momenetum of middle particle - - final_time = num_steps * dt - t_span = (0.0, final_time) - prob = ODEProblem(tendency!, s_0, t_span, params) - sol = solve(prob, SSPRK33(), dt = dt) # perturbed simulation - state_perturbed = sol.u - time_perturbed = sol.t - - # println(time_perturbed .- timesteps, "\n") - - println("Finished simulations") - println("Control len: ", length(states_control), " Perturbed len: ", length(state_perturbed), " Time len: ", length(time_perturbed)) - - - # # print the difference in the two initial states - # println("Initial State Difference (pertrubed - control): ", state_perturbed[1] .- states_control[1]) - # println("Final State Difference (pertrubed - control): ", state_perturbed[end] .- states_control[end]) - - # plot the difference in the two states over time - differences = [] - for i in 1:length(state_perturbed) - # calculate mode 1 difference - a_k = 0 - amp = sqrt(2 / (N + 1)) - for j in 1:N - diff = abs(state_perturbed[i][j*2-1] - states_control[i][j*2-1]) - a_k += diff * sin((1 * (j - 1) * π) / (N + 1)) - end - a_k *= amp - push!(differences, a_k) - end - plot(time_perturbed, differences, label = "Difference in States", xlabel = "Time", ylabel = "Difference", title = "Absolute Difference in States Over Time", lw = 1.5, color = :blue) - savefig("t/frames/difference-in-states-beta-$beta.png") - # perform a linear regression on the log of the differences - # cut the data set into everythign after t=500 - # find the index where t = 500 - # index = 1 - # for i in 1:length(time_perturbed) - # if time_perturbed[i] >= 500 - # index = i - # break - # end - # end - # differences = differences[index:end] - # time_perturbed = time_perturbed[index:end] - ln_differeces = log.(differences) - (a, b) = linear_regression(time_perturbed, ln_differeces) - a = round(a, digits = 4) - b = round(b, digits = 4) - plot(time_perturbed, differences, label = "ln(difference in mode1)", xlabel = "Time", ylabel = "Difference", title = "Absolute Difference in States Over Time", msw = 0.0, yscale = :ln, color = :blue, lw = 1.5) - plot!(time_perturbed, exp.(a * time_perturbed .+ b), label = "exp($a*t + $b)", lw = 2, color = :red, linestyle = :dash) - savefig("t/frames/difference-in-states-beta-$beta-log.png") - - println("Saved difference-in-states-beta-$beta.png and difference-in-states-beta-$beta-log.png\n") - println("Lynapov Exponent: ", a) - - # check that solve doesn't explode - # animate_states(state_perturbed, time_perturbed, 0, 1) - return a -end - -# linear regression of the log of the difference -function linear_regression(x, y) - n = length(x) - x̄ = sum(x) / n - ȳ = sum(y) / n - a = sum((x .- x̄) .* (y .- ȳ)) / sum((x .- x̄) .^ 2) - b = ȳ - a * x̄ - return (a, b) -end - -function analyze_equiparition(t_max, beta) - # run a simulation with normal parameters - N = 32 # number of masses - A = 10 # amplitude - modes = 3 # number of modes to plot - final_time = t_max # seconds - dt = 0.05 # seconds - num_steps = Int(final_time / dt) - - t_span = (0.0, final_time) - params = (N, modes, beta, A, dt, num_steps) - s = run_simulation(N, modes, beta, A, dt, num_steps) - - states = s.u - timesteps = s.t - - analyze_energies_of_n_modes([1, 3, 5, 7], states, timesteps, beta) -end - -# my_beta = 0.3 -# estimate_lynapov_exponent(10000, my_beta, 0.001) -# analyze_equiparition(10000, my_beta) - -# function plot_betas_versus_lynapov_exponent(betas, t_max, jitter = 0.001) -# lynapovs = [] -# for b in betas -# lynapov = estimate_lynapov_exponent(t_max, b, jitter) -# push!(lynapovs, lynapov) -# end - -# plot(betas, lynapovs, label = "Lynapov Exponent", xlabel = "Beta", ylabel = "Lynapov Exponent", title = "Lynapov Exponent Over Beta") -# savefig("t/lynapov-over-beta.png") -# end - -function plot_probability_distrubtion_of_mode(energy_of_mode, total_energy, timesteps, beta) - # get the probability distribution of the mode - prob = [] - for i in 1:length(timesteps) - push!(prob, energy_of_mode[i] / total_energy[1]) # energy is conserved - end - - plot(timesteps, prob, label = "Probability of Mode 1", xlabel = "Time", ylabel = "Probability", title = "Probability of Mode 1 Over Time with beta=$beta", ylim = (0, 1)) - savefig("t/probability-mode1-beta-$beta.png") -end - -function analyze_probability_distribution_of_mode(t_max, beta) - # run a simulation with normal parameters - N = 32 # number of masses - A = 10 # amplitude - modes = 3 # number of modes to plot - final_time = t_max # seconds - dt = 0.05 # seconds - num_steps = Int(final_time / dt) - - t_span = (0.0, final_time) - params = (N, modes, beta, A, dt, num_steps) - s = run_simulation(N, modes, beta, A, dt, num_steps) - - states = s.u - timesteps = s.t - - # get the energies of the mode - total_energy, kinetic_energy, potential_energy, total_by_mode, kinetic_by_mode, potential_by_mode = analyze_energies_of_n_modes([1, 3, 5, 7, 9, 11, 13], states, timesteps, beta) - - # get the probability distribution of the mode - # plot_probability_distrubtion_of_mode(total_by_mode[1], total_energy, timesteps, beta) - - # print the time when the probability gets lower than .5 - for i in 1:length(timesteps) - if (total_by_mode[1][i] / total_energy[1]) < 0.6 - println("TIME when probability of mode 1 is less than 0.6 for beta=$beta: ", timesteps[i]) - return timesteps[i] - end - end - - return -1 -end - -# estimate_lynapov_exponent(1000, 0.3, 0.001) - -function plot_pos_at_time_t(t, beta) - # run a simulation with normal parameters - N = 32 # number of masses - A = 10 # amplitude - modes = 3 # number of modes to plot - final_time = t # seconds - dt = 0.05 # seconds - num_steps = Int(final_time / dt) - - t_span = (0.0, final_time) - params = (N, modes, beta, A, dt, num_steps) - s = run_simulation(N, modes, beta, A, dt, num_steps) - - states = s.u - timesteps = s.t - positions = s.u[end][1:2:end] - inital_positions = s.u[1][1:2:end] - - plot(positions, label = "Positions at t=$t", xlabel = "Mass Number", ylabel = "Displacement", title = "Positions of Masses at t=$t, beta=$beta", lw = 2, marker = :circle) - plot!(inital_positions, label = "Initial Positions", lw = 2, marker = :circle) - savefig("t/new/positions-at-t-$t-beta-$beta.png") -end - -function get_breakdown_time(beta, t_find = 3600) - t = analyze_probability_distribution_of_mode(t_find, beta) - if t == -1 - println("Could not find time for beta: ", beta) - return -1 - end - return t -end - - -# plot_pos_at_time_t(1000, 10) -# plot_pos_at_time_t(1000, 0) - - - -# plot_betas_versus_lynapov_exponent(collect(0:1:50), 1000, 0.001) - -# analyze_probability_distribution_of_mode(100000, 2.0) - -# break_down_times = [] -# bs = collect(1:0.5:50) -# for b_s in bs -# t_find = 3600 -# t = analyze_probability_distribution_of_mode(t_find, b_s) -# if t == -1 -# println("Could not find time for beta: ", b_s) -# continue -# end -# push!(break_down_times, t) -# end -# # plot beta vs breakdown time -# plot(bs, break_down_times, label = "Breakdown Time", xlabel = "Beta", ylabel = "Time", title = "Breakdown Time Over Beta") -# savefig("t/breakdown-time-over-beta.png") - -# # plot the log of this -# one_over_breakdown_time = break_down_times .^ (-1) -# (a, b) = linear_regression(bs, one_over_breakdown_time) -# a = round(a, digits = 4) -# b = round(b, digits = 4) -# plot(bs, one_over_breakdown_time, label = "1 / (breakdown time)", xlabel = "Beta", ylabel = "1 / (breakdown time)", title = "Breakdown Time ^ -1 Over Beta", msw = 0.0, color = :blue, lw = 1.5) -# plot!(bs, a * bs .+ b, label = "1 / t = $a * beta + $b", lw = 2, color = :red, linestyle = :dash) -# savefig("t/ln-breakdown-time-over-beta.png") - -# println("Done!") -# # for a range of betas, find the breakdown time and the lynapov exponent -function tmp(betas) - break_down_times = [] - lynapovs = [] - t_find = 350000 - for b in betas - # if b < 0.5 - # t_find = 10000 - # end - t = get_breakdown_time(b, t_find) - if t == -1 - continue - end - t_find = round(Int, t) + 100 - - l = estimate_lynapov_exponent(t_find, b, 0.001) - if l == 0 - println("Lynapov for beta was zero: ", b) - continue - end - push!(break_down_times, t) - push!(lynapovs, l) - end - - # plot lynapov vs breakdown time - scatter(break_down_times, lynapovs, label = "Lynapov Exponent", xlabel = "Breakdown Time", ylabel = "Lynapov Exponent", title = "Lynapov Exponent Over Breakdown Time") - savefig("t/new/new-lynapov-over-breakdown-time.png") - - # plot this over log of lynapovs - ln_lynapovs = log.(lynapovs) - ln_break_down_times = log.(break_down_times) - (a, b) = linear_regression(ln_break_down_times, ln_lynapovs) - a = round(a, digits = 4) - b = round(b, digits = 4) - scatter(ln_break_down_times, ln_lynapovs, label = "ln(Lynapov Exponent)", xlabel = "ln(Breakdown Time)", ylabel = "ln(Lynapov Exponent)", title = "ln(Lynapov Exponent) Over ln(Breakdown Time)", msw = 0.0, color = :blue, lw = 1.5) - plot!(ln_break_down_times, a * ln_break_down_times .+ b, label = "ln(beta) = $a * ln(t) + $b", lw = 2, color = :red, linestyle = :dash) - savefig("t/new/new-ln-lynapov-over-breakdown-time.png") - - # only one log - (a, b) = linear_regression(lynapovs, ln_break_down_times) - a = round(a, digits = 4) - b = round(b, digits = 4) - scatter(lynapovs, ln_break_down_times, label = "ln(Breakdown Time)", xlabel = "Lynapov Exponent", ylabel = "ln(Breakdown Time)", title = "ln(Breakdown Time) Over Lynapov Exponent", msw = 0.0, color = :blue, lw = 1.5) - plot!(lynapovs, a * lynapovs .+ b, label = "ln(t) = $a * beta + $b", lw = 2, color = :red, linestyle = :dash) - savefig("t/new/new-new-ln-breakdown-time-over-lynapov.png") - - - - # plot the lynapoovs to the ^-1 - one_over_breakdown_time = break_down_times .^ (-1) - (a, b) = linear_regression(lynapovs, one_over_breakdown_time) - a = round(a, digits = 4) - b = round(b, digits = 4) - scatter(lynapovs, one_over_breakdown_time, label = "1 / (breakdown time)", xlabel = "Lynapov Exponent", ylabel = "1 / (breakdown time)", title = "1 / (breakdown time) Over Lynapov Exponent", msw = 0.0, color = :blue, lw = 1.5) - plot!(lynapovs, a * lynapovs .+ b, label = "1 / t = $a * beta + $b", lw = 2, color = :red, linestyle = :dash) - savefig("t/new/new-one-over-lynapov-over-breakdown-time.png") - - # plot the lynapov over betas - scatter(betas, lynapovs, label = "Lynapov Exponent", xlabel = "Beta", ylabel = "Lynapov Exponent", title = "Lynapov Exponent Over Beta") - savefig("t/new/new-lynapov-over-beta.png") -end - -tmp(collect(0.3:0.025:10)) - -# estimate_lynapov_exponent(500000, 0.3, 0.001) |