diff options
Diffstat (limited to 'hw3/Lorenz.jl')
| -rw-r--r-- | hw3/Lorenz.jl | 80 |
1 files changed, 44 insertions, 36 deletions
diff --git a/hw3/Lorenz.jl b/hw3/Lorenz.jl index c4c501d..4a7fa0e 100644 --- a/hw3/Lorenz.jl +++ b/hw3/Lorenz.jl @@ -2,7 +2,7 @@ using DifferentialEquations using Plots r_max = 160 -r_steps = 320 +dr = .01 sim_time = 1000.0 @@ -19,49 +19,57 @@ end # make a linspace for these values of r +r_steps = Int(r_max/dr) r_values = range(0, r_max, length=r_steps) -sols = [] -for i in 1:r_steps - r = r_values[i] - p = [10.0, r, 8/3] # parameters of the Lorentz 63 system - tspan = (0.0, sim_time) - - u0 = [1.0, 0.0, 0.0] - prob = ODEProblem(tendency!, u0, tspan, p) - sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # control simulation - - push!(sols, sol) - - println("r=$r") -end - -# STEP 2, map data to arrays where plane crosses the x-axis r_maxes = [] z_maxes = [] r_mins = [] z_mins = [] -for i in 1:r_steps - println("i: ", length(sols[i].t)) - z_values = sols[i][3, :] - # iterate over to find the local maxima and minima - # take off the first 300 values to avoid transient - for j in 301:length(z_values)-1 - if z_values[j] > z_values[j-1] && z_values[j] > z_values[j+1] - push!(r_maxes, r_values[i]) - push!(z_maxes, z_values[j]) +let + r_diverge_1 = -1 + for i in 1:r_steps + r = r_values[i] + p = [10.0, r, 8/3] # parameters of the Lorentz 63 system + tspan = (0.0, sim_time) + + u0 = [1.0, 0.0, 0.0] + prob = ODEProblem(tendency!, u0, tspan, p) + sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # control simulation + + # iterate over to find the local maxima and minima + z_values = sol[3, :] + + # take off the first 300 values to avoid transient + for j in 301:length(z_values)-1 + if z_values[j] > z_values[j-1] && z_values[j] > z_values[j+1] + push!(r_maxes, r_values[i]) + push!(z_maxes, z_values[j]) + end + if z_values[j] < z_values[j-1] && z_values[j] < z_values[j+1] + push!(r_mins, r_values[i]) + push!(z_mins, z_values[j]) + end + + # calcualte dz + # if r_diverge_1 < 0 && length(z_mins) > 1 && length(z_maxes) > 1 && abs(z_mins[end] - z_maxes[end]) > 20.0 + # r_diverge_1 = r + # println("Divergence starts at r = $r") + # end end - if z_values[j] < z_values[j-1] && z_values[j] < z_values[j+1] - push!(r_mins, r_values[i]) - push!(z_mins, z_values[j]) + + if i % 500 == 0 + println("r=$r") end end -end -# println("r_maxes: ", r_maxes) -# println("z_maxes: ", z_maxes) + plt = plot(r_maxes, z_maxes, seriestype=:scatter, mc=:blue, ms=.25, ma=.25, msw=0, xlabel="r", ylabel="z", title="Bifurcation Diagram for Lorenz 63 System", legend=:topleft, grid=:true, label="Local Maxima") + plot!(r_mins, z_mins, seriestype=:scatter, mc=:green, ms=.25, ma=.25, msw=0, xlabel="r", ylabel="z", label="Local Minima") + + # r_diverge_display = round(r_diverge_1, digits=3) + vline!([24], label="Critical r value ≈ 24", lc=:red, lw=1.5, ls=:dash) + vline!([122], label="Critical r value ≈ 122", lc=:orange, lw=1.5, ls=:dash) + vline!([146], label="Critical r value ≈ 146", lc=:purple, lw=1.5, ls=:dash) -# STEP 3, plot the bifurcation diagram -plot(r_maxes, z_maxes, seriestype = :scatter, mc=:blue, ms=.25, ma=0.25, label="Maxima") -plot!(r_mins, z_mins, seriestype = :scatter, mc=:green, ms=.25, ma=0.25, label="Minima") -savefig("hw3/test3.png")
\ No newline at end of file + savefig("hw3/test11.png") +end
\ No newline at end of file |