From e650ed1e1e908e51c78c1b047bec0da7c4fea366 Mon Sep 17 00:00:00 2001 From: sotech117 Date: Sat, 27 Apr 2024 04:25:23 -0400 Subject: testing --- hw7/.DS_Store | Bin 0 -> 8196 bytes hw7/5.6-face.png | Bin 0 -> 76602 bytes hw7/5.6-fx.png | Bin 0 -> 38040 bytes hw7/5.6-fy.png | Bin 0 -> 64046 bytes hw7/5.6-phi-field.png | Bin 0 -> 25715 bytes hw7/5.6-phi.png | Bin 0 -> 54394 bytes hw7/5.6-test.png | Bin 0 -> 26702 bytes hw7/5.6.jl | 192 +++++++++++++++++++++++++++++++++ hw7/5.6.png | Bin 0 -> 86224 bytes hw7/8-12.jl | 191 +++++++++++++++++++++++++++++++++ hw7/8-15-d.png | Bin 0 -> 30807 bytes hw7/8-15-scales-d.png | Bin 0 -> 28596 bytes hw7/8-15-scales.png | Bin 0 -> 31626 bytes hw7/8-15.jl | 274 +++++++++++++++++++++++++++++++++++++++++++++++ hw7/8-15.png | Bin 0 -> 31061 bytes hw7/BoundStates.jl | 73 +++++++++++++ hw7/Laplacians.jl | 145 +++++++++++++++++++++++++ hw7/RadialBoundStates.jl | 59 ++++++++++ hw7/mvsb.svg | 55 ++++++++++ hw7/mvsb2.svg | 46 ++++++++ hw7/tightbinding (1).jl | 24 +++++ 21 files changed, 1059 insertions(+) create mode 100644 hw7/.DS_Store create mode 100644 hw7/5.6-face.png create mode 100644 hw7/5.6-fx.png create mode 100644 hw7/5.6-fy.png create mode 100644 hw7/5.6-phi-field.png create mode 100644 hw7/5.6-phi.png create mode 100644 hw7/5.6-test.png create mode 100644 hw7/5.6.jl create mode 100644 hw7/5.6.png create mode 100644 hw7/8-12.jl create mode 100644 hw7/8-15-d.png create mode 100644 hw7/8-15-scales-d.png create mode 100644 hw7/8-15-scales.png create mode 100644 hw7/8-15.jl create mode 100644 hw7/8-15.png create mode 100644 hw7/BoundStates.jl create mode 100644 hw7/Laplacians.jl create mode 100644 hw7/RadialBoundStates.jl create mode 100644 hw7/mvsb.svg create mode 100644 hw7/mvsb2.svg create mode 100644 hw7/tightbinding (1).jl (limited to 'hw7') diff --git a/hw7/.DS_Store b/hw7/.DS_Store new file mode 100644 index 0000000..9b9f923 Binary files /dev/null and b/hw7/.DS_Store differ diff --git a/hw7/5.6-face.png b/hw7/5.6-face.png new file mode 100644 index 0000000..b02d8f9 Binary files /dev/null and b/hw7/5.6-face.png differ diff --git a/hw7/5.6-fx.png b/hw7/5.6-fx.png new file mode 100644 index 0000000..a445d09 Binary files /dev/null and b/hw7/5.6-fx.png differ diff --git a/hw7/5.6-fy.png b/hw7/5.6-fy.png new file mode 100644 index 0000000..56c38f5 Binary files /dev/null and b/hw7/5.6-fy.png differ diff --git a/hw7/5.6-phi-field.png b/hw7/5.6-phi-field.png new file mode 100644 index 0000000..3556e4b Binary files /dev/null and b/hw7/5.6-phi-field.png differ diff --git a/hw7/5.6-phi.png b/hw7/5.6-phi.png new file mode 100644 index 0000000..93b9775 Binary files /dev/null and b/hw7/5.6-phi.png differ diff --git a/hw7/5.6-test.png b/hw7/5.6-test.png new file mode 100644 index 0000000..58328e6 Binary files /dev/null and b/hw7/5.6-test.png differ diff --git a/hw7/5.6.jl b/hw7/5.6.jl new file mode 100644 index 0000000..e6f6c8a --- /dev/null +++ b/hw7/5.6.jl @@ -0,0 +1,192 @@ +using Pkg +Pkg.add("Plots") +using Plots + +function del2_5(a::Matrix{Float64}, dx = 1.0) + #= + Returns a finite-difference approximation of the laplacian of the array a, + with lattice spacing dx, using the five-point stencil: + + 0 1 0 + 1 -4 1 + 0 1 0 + =# + + del2 = zeros(size(a)) + del2[2:end-1, 2:end-1] .= (a[2:end-1, 3:end] + a[2:end-1, 1:end-2] + + a[3:end, 2:end-1] + a[1:end-2, 2:end-1] - + 4.0 * a[2:end-1, 2:end-1]) ./ (dx^2) + return del2 +end + +function del2_9(a::Matrix{Float64}, dx = 1.0) + #= + Returns a finite-difference approximation of the laplacian of the array a, + with lattice spacing dx, using the nine-point stencil: + + 1/6 2/3 1/6 + 2/3 -10/3 2/3 + 1/6 2/3 1/6 + =# + + del2 = zeros(size(a)) + del2[2:end-1, 2:end-1] .= (4.0 * (a[2:end-1, 3:end] + a[2:end-1, 1:end-2] + + a[3:end, 2:end-1] + a[1:end-2, 2:end-1]) + + (a[1:end-2, 1:end-2] + a[1:end-2, 3:end] + + a[3:end, 1:end-2] + a[3:end, 3:end]) - + 20.0 * a[2:end-1, 2:end-1]) ./ (6.0 * dx^2) + return del2 +end + +function invDel2_5(b::Matrix{Float64}, dx = 1.0, N = 10000) + #= + Relaxes over N steps to a discrete approximation of the inverse laplacian + of the source term b, each step is a weighted average over the four neighboring + points and the source. This is the Jacobi algorithm. + =# + + invDel2 = zeros(size(b)) + newInvDel2 = zeros(size(b)) + + for m in 1:N + newInvDel2[2:end-1, 2:end-1] .= 0.25 * (invDel2[2:end-1, 3:end] + + invDel2[2:end-1, 1:end-2] + + invDel2[3:end, 2:end-1] + + invDel2[1:end-2, 2:end-1] - + (dx^2) * b[2:end-1, 2:end-1]) + invDel2 .= newInvDel2 + end + + diff = del2_5(invDel2, dx) - b + diffSq = diff .* diff + error = sqrt(sum(diffSq)) + + println("\nerror = ", error) + + return invDel2 +end + + +# # Define variables +# L = 10 +# dx = 1.0 + +# # Parabola of revolution has constant Laplacian +# phi = zeros((L, L)) +# for i ∈ 1:L +# x = (i + 0.5 - 0.5 * L) * dx +# for j ∈ 1:L +# y = (j + 0.5 - 0.5 * L) * dx +# phi[i, j] = x^2 + y^2 +# end +# end + +# println(size(phi)) + +# println(del2_5(phi, dx)) + +# println("\n\n") + +# println(del2_9(phi, dx)) + +# Electrostatics example: uniformly charged cylinder (rho = 1) of radius R +L = 100 +dx = 1.0 +R = 1.0 +R2 = R^2 + +# Charge distribution +rho = zeros((L, L)) +for i ∈ 1:L + y = (i - 1 + 0.5 - 0.5 * L) * dx + println("y = ", y) + for j ∈ 1:L + x = (j - 1 + 0.5 - 0.5 * L) * dx + # this is the rod here + if y < 0 && abs(x) < R + rho[i, j] = 1.0 # high voltage + elseif y > 0 && abs(x) < R + rho[i, j] = -1.0 # conducting plane + else + rho[i, j] = 0.0 + end + end +end + +rhoPlot = plot(rho[Int(L / 2), :]) + +# Exact (analytical) electric potential +# phi = zeros((L, L)) +# for i ∈ 1:L +# x = (i + 0.5 - 0.5 * L) * dx +# for j ∈ 1:L +# y = (j + 0.5 - 0.5 * L) * dx +# r2 = x^2 + y^2 +# if r2 < R2 +# phi[i, j] = -π * r2 +# else +# phi[i, j] = -π * (R2 + R2 * log(r2 / R2)) +# end +# end +# end + +# phiPlot = plot(phi[Int(L / 2), :]) + +# Charge density obtained from exact potential +# rho = -1.0 / (4.0 * π) * del2_5(phi, dx) +# rhoPlotLattice = plot(rho[Int(L / 2), :]) + +function calc_field(phi::Matrix{Float64}, dx = 1.0) + E_x = zeros(size(phi)) + E_y = zeros(size(phi)) + + E_x[2:end-1, 2:end-1] .= -(phi[3:end, 2:end-1] - phi[1:end-2, 2:end-1]) / (2 * dx) + E_y[2:end-1, 2:end-1] .= -(phi[2:end-1, 3:end] - phi[2:end-1, 1:end-2]) / (2 * dx) + + # for i ∈ 2:size(phi, 1)-1 + # for j ∈ 2:size(phi, 2)-1 + # E_x[i, j] = -(phi[i+1, j] - phi[i-1, j]) / (2 * dx) + # E_y[i, j] = -(phi[i, j+1] - phi[i, j-1]) / (2 * dx) + # end + # end + + return E_x, E_y +end + +phi = invDel2_5(-4.0 * π * rho, dx, 1000) +# phi = invDelSOR(-4.0 * π * rho, dx, 500) +phi .-= phi[Int(L / 2), Int(L / 2)] + +phiPlotInvDel = plot(phi[Int(L / 2), :]) + +# phi_contour = contourf(phi, title = "Electric Potential", color = :viridis, aspect_ratio = :equal, colorbar_title = "V") +# savefig(phi_contour, "hw7/5.6-phi.png") + +# plot the arrows of the magnatic field on the contour plot +xxs = [x for x in 1:L for y in 1:L] +xxy = [y for x in 1:L for y in 1:L] +d_x, d_y = calc_field(phi, dx) + +function df(i, j) + norm = d_x[Int64(i)] * d_x[Int64(i)] + d_y[Int64(j)] * d_y[Int64(j)] + if norm == 0 + return (0, 0) + end + return (d_x[Int64(i)] / norm, d_y[Int64(j)] / norm) +end + +quiver(xxs, xxy, quiver = df, aspect_ratio = :equal, legend = false) +savefig("hw7/5.6-test.png") + + +# d_x, d_y = calc_field(phi, dx) +# field_y = contourf(d_y, title = "Electric Field (y)", color = :viridis, aspect_ratio = :equal, colorbar_title = "V/m") + +# field_x = contourf(d_x, title = "Electric Field (x)", color = :viridis, aspect_ratio = :equal, colorbar_title = "V/m") + + +# savefig(field_y, "hw7/5.6-fy.png") +# savefig(field_x, "hw7/5.6-fx.png") + +#plot(rhoPlot, rhoPlotLattice) +#plot(phiPlot, phiPlotInvDel) diff --git a/hw7/5.6.png b/hw7/5.6.png new file mode 100644 index 0000000..f17824a Binary files /dev/null and b/hw7/5.6.png differ diff --git a/hw7/8-12.jl b/hw7/8-12.jl new file mode 100644 index 0000000..29e18ec --- /dev/null +++ b/hw7/8-12.jl @@ -0,0 +1,191 @@ +#!/Applications/Julia-1.7.app/Contents/Resources/julia/bin/julia + +using Statistics +using Plots + +function wrap_index(i::Int, l::Int)::Int + wrap = (i - 1) % l + 1 + return (wrap <= 0) ? l + wrap : wrap +end + +mutable struct Ising2D + l::Int + n::Int + temperature::Float64 + w::Vector{Float64} # Boltzmann weights + state::Matrix + energy::Float64 + magnetization::Int + mc_steps::Int + accepted_moves::Int + energy_array::Vector{Float64} + magnetization_array::Vector{Int} + H::Float64 +end + +Ising2D(l::Int, temperature::Float64, H=1.0) = begin + n = l^2 + w = zeros(9) + w[9] = exp(-8.0 / temperature) + w[5] = exp(-4.0 / temperature) + state = ones(Int, l, l) # initially all spins up + energy = Float64(-2 * n + 2 * H * n) + magnetization = n + return Ising2D(l, n, temperature, w, state, energy, magnetization, 0, 0, + Int[], Int[], H) +end + +function reset!(ising::Ising2D) + ising.mc_steps = 0 + ising.accepted_moves = 0 + ising.energy_array = Int[] + ising.magnetization_array = Int[] +end + +function mc_step!(ising::Ising2D) + l::Int = ising.l + n::Int = ising.n + w = ising.w + + state = ising.state + accepted_moves = ising.accepted_moves + energy = ising.energy + magnetization = ising.magnetization + + random_positions = l * rand(2 * n) + random_array = rand(n) + + for k in 1:n + i = trunc(Int, random_positions[2 * k - 1]) + 1 + j = trunc(Int, random_positions[2 * k]) + 1 + + changed_spins = state[i, j] * (state[i % l + 1, j] + + state[wrap_index(i - 1, l), j] + state[i, j % l + 1] + + state[i, wrap_index(j - 1, l)]) + de = 2 * changed_spins + 2 * ising.H * state[i, j] + + if de <= 0 || rand() < exp(-de / ising.temperature) + accepted_moves += 1 + new_spin = - state[i, j] # flip spin + state[i, j] = new_spin + + # add the effects of the new spin + energy += de + magnetization += 2 * new_spin + end + + end + + ising.state = state + ising.accepted_moves = accepted_moves + ising.energy = energy + ising.magnetization = magnetization + + append!(ising.energy_array, ising.energy) + append!(ising.magnetization_array, ising.magnetization) + ising.mc_steps = ising.mc_steps + 1 +end + +function steps!(ising::Ising2D, num::Int=100) + for i in 1:num + mc_step!(ising) + end +end + +function mean_energy(ising::Ising2D) + return mean(ising.energy_array) / ising.n +end + +function specific_heat(ising::Ising2D) + return (std(ising.energy_array) / ising.temperature) ^ 2 / ising.n +end + +function mean_magnetization(ising::Ising2D) + return mean(ising.magnetization_array) / ising.n +end + +function susceptibility(ising::Ising2D) + return (std(ising.magnetization_array)) ^ 2 / (ising.temperature * ising.n) +end + +function observables(ising::Ising2D) + printstyled("Temperature\t\t", bold=true) + print(ising.temperature); print("\n") + + printstyled("Mean Energy\t\t", bold=true) + print(mean_energy(ising)); print("\n") + + printstyled("Mean Magnetiz.\t\t", bold=true) + print(mean_magnetization(ising)); print("\n") + + printstyled("Specific Heat\t\t", bold=true) + print(specific_heat(ising)); print("\n") + + printstyled("Susceptibility\t\t", bold=true) + print(susceptibility(ising)); print("\n") + + printstyled("MC Steps\t\t", bold=true) + print(ising.mc_steps); print("\n") + printstyled("Accepted Moves\t\t", bold=true) + print(ising.accepted_moves); print("\n") +end + + +function plot_ising(state::Matrix{Int}) + pos = Tuple.(findall(>(0), state)) + neg = Tuple.(findall(<(0), state)) + scatter(pos, markersize=5) + scatter!(neg, markersize=5) +end + +function find_m(H::Float64, l::Int, num::Int, T::Float64) + m = Ising2D(l, T, H) + steps!(m, num) + print("T = $T\n") + print("H = $H\n") + print("Mean Energy: $(mean_energy(m))\n") + print("Mean Magnetization: $(mean_magnetization(m))\n\n") + return mean_magnetization(m) +end + +function map_h_to_m(H_range::Vector{Float64}, l::Int, num::Int, T::Float64) + m = [] + for H in H_range + push!(m, find_m(H, l, num, T)) + end + return m +end + +function do_linear_regression(x::Vector{Float64}, y::Vector{Float64}) + n = length(x) + x̄ = mean(x) + ȳ = mean(y) + Σxy = sum((x .- x̄) .* (y .- ȳ)) + Σx² = sum((x .- x̄) .^ 2) + b = Σxy / Σx² + a = ȳ - b * x̄ + return a, b +end + +function plot_log_of_m_and_h(H_range::Vector{Float64}, l::Int, num::Int, T=2.27) + m = map_h_to_m(H_range, l, num, T) + p = scatter(H_range, m, label="M vs H", xlabel="H", ylabel="M", title="Magnetization (M) vs Field (B) for Ising Model at T_c", scale=:ln) + + # get the linear regression of the log + log_h = log.(H_range) + log_m = log.(m) + a, b = do_linear_regression(log_h, log_m) + println("a: $a, b: $b") + # plot the linear regression + plot!(p, H_range, exp.(a) .* H_range .^ b, label="linear regression = $(round(a, digits=3)) + $(round(b, digits=3))x", line=:dash, color=:red) + + return p +end + +# textbook rec +h_range = .02:.02:.2 +h_range = collect(h_range) +T = 2.27 # T_c for this system +side = 64 +steps = 5000 +plot_log_of_m_and_h(h_range, side, steps) \ No newline at end of file diff --git a/hw7/8-15-d.png b/hw7/8-15-d.png new file mode 100644 index 0000000..eb64f61 Binary files /dev/null and b/hw7/8-15-d.png differ diff --git a/hw7/8-15-scales-d.png b/hw7/8-15-scales-d.png new file mode 100644 index 0000000..5e4f2e4 Binary files /dev/null and b/hw7/8-15-scales-d.png differ diff --git a/hw7/8-15-scales.png b/hw7/8-15-scales.png new file mode 100644 index 0000000..2256167 Binary files /dev/null and b/hw7/8-15-scales.png differ diff --git a/hw7/8-15.jl b/hw7/8-15.jl new file mode 100644 index 0000000..5c833be --- /dev/null +++ b/hw7/8-15.jl @@ -0,0 +1,274 @@ +#!/Applications/Julia-1.7.app/Contents/Resources/julia/bin/julia + +using Statistics +using Plots + +function wrap_index(i::Int, l::Int)::Int + wrap = (i - 1) % l + 1 + return (wrap <= 0) ? l + wrap : wrap +end + +mutable struct Ising2D + l::Int + n::Int + temperature::Float64 + w::Vector{Float64} # Boltzmann weights + state::Matrix + energy::Float64 + magnetization::Int + mc_steps::Int + accepted_moves::Int + energy_array::Vector{Float64} + magnetization_array::Vector{Int} + H::Float64 +end + +Ising2D(l::Int, temperature::Float64, H = 1.0) = begin + n = l^2 + w = zeros(9) + w[9] = exp(-8.0 / temperature) + w[5] = exp(-4.0 / temperature) + state = ones(Int, l, l) # initially all spins up + energy = Float64(-2 * n + 2 * H * n) + magnetization = n + return Ising2D(l, n, temperature, w, state, energy, magnetization, 0, 0, + Int[], Int[], H) +end + +function reset!(ising::Ising2D) + ising.mc_steps = 0 + ising.accepted_moves = 0 + ising.energy_array = Int[] + ising.magnetization_array = Int[] +end + +function mc_step!(ising::Ising2D) + l::Int = ising.l + n::Int = ising.n + w = ising.w + + state = ising.state + accepted_moves = ising.accepted_moves + energy = ising.energy + magnetization = ising.magnetization + + random_positions = l * rand(2 * n) + random_array = rand(n) + + for k in 1:n + i = trunc(Int, random_positions[2*k-1]) + 1 + j = trunc(Int, random_positions[2*k]) + 1 + + changed_spins = state[i, j] * (state[i%l+1, j] + + state[wrap_index(i - 1, l), j] + state[i, j%l+1] + + state[i, wrap_index(j - 1, l)]) + de = 2 * changed_spins + 2 * ising.H * state[i, j] + + if de <= 0 || rand() < exp(-de / ising.temperature) + accepted_moves += 1 + new_spin = -state[i, j] # flip spin + state[i, j] = new_spin + + # add the effects of the new spin + energy += de + magnetization += 2 * new_spin + end + + end + + ising.state = state + ising.accepted_moves = accepted_moves + ising.energy = energy + ising.magnetization = magnetization + + append!(ising.energy_array, ising.energy) + append!(ising.magnetization_array, ising.magnetization) + ising.mc_steps = ising.mc_steps + 1 +end + +function steps!(ising::Ising2D, num::Int = 100) + for i in 1:num + mc_step!(ising) + end +end + +function mean_energy(ising::Ising2D) + return mean(ising.energy_array) / ising.n +end + +function specific_heat(ising::Ising2D) + return (std(ising.energy_array) / ising.temperature)^2 / ising.n +end + +function mean_magnetization(ising::Ising2D) + return mean(ising.magnetization_array) / ising.n +end + +function susceptibility(ising::Ising2D) + return (std(ising.magnetization_array))^2 / (ising.temperature * ising.n) +end + +function observables(ising::Ising2D) + printstyled("Temperature\t\t", bold = true) + print(ising.temperature) + print("\n") + + printstyled("Mean Energy\t\t", bold = true) + print(mean_energy(ising)) + print("\n") + + printstyled("Mean Magnetiz.\t\t", bold = true) + print(mean_magnetization(ising)) + print("\n") + + printstyled("Specific Heat\t\t", bold = true) + print(specific_heat(ising)) + print("\n") + + printstyled("Susceptibility\t\t", bold = true) + print(susceptibility(ising)) + print("\n") + + printstyled("MC Steps\t\t", bold = true) + print(ising.mc_steps) + print("\n") + printstyled("Accepted Moves\t\t", bold = true) + print(ising.accepted_moves) + print("\n") +end + + +function plot_ising(state::Matrix{Int}) + pos = Tuple.(findall(>(0), state)) + neg = Tuple.(findall(<(0), state)) + scatter(pos, markersize = 5) + scatter!(neg, markersize = 5) +end + +function h_to_m(H::Float64, l::Int, num::Int, T::Float64) + m = Ising2D(l, T, H) + steps!(m, num) + print("T = $T\n") + print("H = $H\n") + print("Mean Energy: $(mean_energy(m))\n") + print("Mean Magnetization: $(mean_magnetization(m))\n\n") + return mean_magnetization(m) +end + +function map_h_to_m(H_range::Vector{Float64}, l::Int, num::Int, T::Float64) + m = [] + for H in H_range + push!(m, h_to_m(H, l, num, T)) + end + return m +end + +function do_linear_regression(x::Vector{Float64}, y::Vector{Float64}) + n = length(x) + x̄ = mean(x) + ȳ = mean(y) + Σxy = sum((x .- x̄) .* (y .- ȳ)) + Σx² = sum((x .- x̄) .^ 2) + b = Σxy / Σx² + a = ȳ - b * x̄ + return a, b +end + +function plot_log_of_m_and_h(H_range::Vector{Float64}, l::Int, num::Int, T = 2.27) + m = map_h_to_m(H_range, l, num, T) + p = scatter(H_range, m, label = "M vs H", xlabel = "H", ylabel = "M", title = "Magnetization (M) vs Field (B) for Ising Model at T_c", scale = :ln) + + # get the linear regression of the log + log_h = log.(H_range) + log_m = log.(m) + a, b = do_linear_regression(log_h, log_m) + println("a: $a, b: $b") + # plot the linear regression + plot!(p, H_range, exp.(a) .* H_range .^ b, label = "linear regression = $(round(a, digits=3)) + $(round(b, digits=3))x", line = :dash, color = :red) + + return p +end + +function t_to_m(T::Float64, l::Int, num::Int, H::Float64) + m = Ising2D(l, T, H) + steps!(m, num) + print("T = $T\n") + print("H = $H\n") + print("Mean Energy: $(mean_energy(m))\n") + print("Mean Sus: $(mean_magnetization(m))\n\n") + return susceptibility(m) +end + +function plot_m_over_t(plt, T_range::Vector{Float64}, l::Int, num::Int, H = 0.0) + m = [] + for T in T_range + push!(m, t_to_m(T, l, num, H)) + end + p = scatter!(plt, T_range, m, label = "H = $H", xlabel = "T", ylabel = "X", title = "Susceptibility (X) vs Temperature (T) on Ising Model") + + return p, m +end + +function plot_m_over_t_and_h(T_range::Vector{Float64}, l::Int, num::Int, H_range::Vector{Float64}) + plt = Plots.scatter() + h = [] + for H in H_range + p, m = plot_m_over_t(plt, T_range, l, num, H) + push!(h, m) + end + + # plot the critical temp as a vertical line + plot!(plt, [2.27, 2.27], [-0.01, 30], label = "T_c = 2.27", line = :dash, color = :red) + + return plt, h +end + +function plot_scales(data, t_range, h_range) + x1 = [] + y1 = [] + x2 = [] + y2 = [] + for i in 1:length(h_range) + h = h_range[i] + for j in 1:length(t_range) + t = t_range[j] + + m = data[i][j] + + scaled_t = abs((t - 2.27) / 2.27) + scaled_m = m * (scaled_t^(7.0 / 4.0)) + scaled_h = h / (scaled_t^(15.0 / 8.0)) + if scaled_h > 30 + continue # dont add if too large + end + + if t > 2.27 + push!(x1, scaled_h) + push!(y1, scaled_m) + else + push!(x2, scaled_h) + push!(y2, scaled_m) + end + end + end + + tmp = scatter(x1, y1, label = "T > T_c", xlabel = "h / abs(t)^(15/8)", ylabel = "X * abs(t)^(7/4)", title = "Susceptibility (X) vs Field (H) for Ising Model") + scatter!(tmp, x2, y2, label = "T < T_c", xlabel = "h / abs(t)^(15/8)", ylabel = "X * abs(t)^(7/4)", title = "Susceptibility (X) vs Field (H) for Ising Model") + return tmp +end + +h_range = 0.01:0.01:0.05 +h_range = collect(h_range) +t_range = 1.5:0.1:3.0 +t_range = collect(t_range) +println("t_range: $t_range") +side = 20 +steps = 3000 +plt, data = plot_m_over_t_and_h(t_range, side, steps, h_range) + +savefig(plt, "hw7/8-15.png") + +p = plot_scales(data, t_range, h_range) +savefig(p, "hw7/8-15-scales.png") + + diff --git a/hw7/8-15.png b/hw7/8-15.png new file mode 100644 index 0000000..a5a4b28 Binary files /dev/null and b/hw7/8-15.png differ diff --git a/hw7/BoundStates.jl b/hw7/BoundStates.jl new file mode 100644 index 0000000..386a5da --- /dev/null +++ b/hw7/BoundStates.jl @@ -0,0 +1,73 @@ +#!/usr/bin/env julia + +"""Find eigenstates and eigenenergies of 1-D quantum problems""" + +using LinearAlgebra +using Plots + +N = 1000 # number of lattice points +L = 20.0 # x runs from -L/2 to L/2 +dx = L / N + +D = zeros(N, N) # discrete laplacian operator +V = zeros(N, N) # potential + +for i in 1:N + D[i, i] = -2.0 +end + +for i in 1:N-1 + D[i, i+1] = 1.0 + D[i+1, i] = 1.0 +end + +#println("\nLattice Laplacian operator") +#println(D) + +function potential(x) + """ The potential energy""" + #return 0.0 # particle in a box + #return 0.5 * x^2 # SHO with the spring constant k = 1 + #return -6.0 * x^2 + 8.0 * x^6 # potential with zero ground state energy + #return 0.1 * x^4 - 2.0 * x^2 + 0.0 * x # double-well potential + return 8*x^6 - 8*x^4 - 4*x^2 + 1 # another double-well potential +end + +for i in 1:N + x = (i + 0.5) * dx - 0.5*L # coordinates of lattice points + V[i, i] = potential(x) +end + +H = -0.5 * dx^(-2.0) * D + V # Hamiltonian. Here m = hbar = 1 + +#println("\nMatrix elements of Hamiltonian = ") +#println(H) + +e, v = eigen(H) # diagonalize Hamiltonian + +println("\nGround state energy = ", e[1]) +println("\n1st excited state energy = ", e[2]) +println("\n2nd excited state energy = ", e[3]) +println("\n3rd excited state energy = ", e[4]) +println("\n4th excited state energy = ", e[5]) + + +gs(x) = exp(-0.5*x^2) # Gaussian that is exact ground state of SHO + + +plot(potential) + + +plot(v[:,1]) +#plot(v[:,2]) +#plot(gs) + +#= +eList = zeros(0) +for i in 1:20 + push!(eList, e[i]) +end + + +bar(eList, orientation = :horizontal) +=# diff --git a/hw7/Laplacians.jl b/hw7/Laplacians.jl new file mode 100644 index 0000000..10f0d48 --- /dev/null +++ b/hw7/Laplacians.jl @@ -0,0 +1,145 @@ +using Pkg +Pkg.add("Plots") +using Plots + +function del2_5(a::Matrix{Float64}, dx=1.0) + #= + Returns a finite-difference approximation of the laplacian of the array a, + with lattice spacing dx, using the five-point stencil: + + 0 1 0 + 1 -4 1 + 0 1 0 + =# + + del2 = zeros(size(a)) + del2[2:end-1, 2:end-1] .= (a[2:end-1, 3:end] + a[2:end-1, 1:end-2] + + a[3:end, 2:end-1] + a[1:end-2, 2:end-1] - + 4.0 * a[2:end-1, 2:end-1]) ./ (dx^2) + return del2 +end + +function del2_9(a::Matrix{Float64}, dx=1.0) + #= + Returns a finite-difference approximation of the laplacian of the array a, + with lattice spacing dx, using the nine-point stencil: + + 1/6 2/3 1/6 + 2/3 -10/3 2/3 + 1/6 2/3 1/6 + =# + + del2 = zeros(size(a)) + del2[2:end-1, 2:end-1] .= (4.0 * (a[2:end-1, 3:end] + a[2:end-1, 1:end-2] + + a[3:end, 2:end-1] + a[1:end-2, 2:end-1]) + + (a[1:end-2, 1:end-2] + a[1:end-2, 3:end] + + a[3:end, 1:end-2] + a[3:end, 3:end]) - + 20.0 * a[2:end-1, 2:end-1]) ./ (6.0 * dx^2) + return del2 +end + +function invDel2_5(b::Matrix{Float64}, dx=1.0, N=10000) + #= + Relaxes over N steps to a discrete approximation of the inverse laplacian + of the source term b, each step is a weighted average over the four neighboring + points and the source. This is the Jacobi algorithm. + =# + + invDel2 = zeros(size(b)) + newInvDel2 = zeros(size(b)) + + for m in 1:N + newInvDel2[2:end-1, 2:end-1] .= 0.25 * (invDel2[2:end-1, 3:end] + + invDel2[2:end-1, 1:end-2] + + invDel2[3:end, 2:end-1] + + invDel2[1:end-2, 2:end-1] - + (dx^2) * b[2:end-1, 2:end-1]) + invDel2 .= newInvDel2 + end + + diff = del2_5(invDel2, dx) - b + diffSq = diff .* diff + error = sqrt(sum(diffSq)) + + println("\nerror = ", error) + + return invDel2 +end + + +# Define variables +L = 10 +dx = 1.0 + +# Parabola of revolution has constant Laplacian +phi = zeros((L, L)) +for i = 1:L + x = (i + 0.5 - 0.5 * L) * dx + for j = 1:L + y = (j + 0.5 - 0.5 * L) * dx + phi[i, j] = x^2 + y^2 + end +end + +println(size(phi)) + +println(del2_5(phi, dx)) + +println("\n\n") + +println(del2_9(phi, dx)) + +# Electrostatics example: uniformly charged cylinder (rho = 1) of radius R +L = 100 +dx = 1.0 +R = 20.0 +R2 = R^2 + +# Charge distribution +rho = zeros((L, L)) +for i = 1:L + x = (i + 0.5 - 0.5 * L) * dx + for j = 1:L + y = (j + 0.5 - 0.5 * L) * dx + r2 = x^2 + y^2 + if r2 < R2 + rho[i, j] = 1.0 + else + rho[i, j] = 0.0 + end + end +end + +rhoPlot = plot(rho[Int(L/2), :]) + +# Exact (analytical) electric potential +phi = zeros((L, L)) +for i = 1:L + x = (i + 0.5 - 0.5 * L) * dx + for j = 1:L + y = (j + 0.5 - 0.5 * L) * dx + r2 = x^2 + y^2 + if r2 < R2 + phi[i, j] = -π * r2 + else + phi[i, j] = -π * (R2 + R2*log(r2/R2)) + end + end +end + +phiPlot = plot(phi[Int(L/2), :]) + +# Charge density obtained from exact potential +rho = -1.0/(4.0 * π) * del2_5(phi, dx) +rhoPlotLattice = plot(rho[Int(L/2), :]) + +phi = invDel2_5(-4.0 * π * rho, dx, 20000) +#phi = invDelSOR(-4.0 * π * rho, dx, 500) +phi .-= phi[Int(L/2), Int(L/2)] + +phiPlotInvDel = plot(phi[Int(L/2), :]) + +contourf(phi) + +#plot(rhoPlot, rhoPlotLattice) +#plot(phiPlot, phiPlotInvDel) diff --git a/hw7/RadialBoundStates.jl b/hw7/RadialBoundStates.jl new file mode 100644 index 0000000..afb91b0 --- /dev/null +++ b/hw7/RadialBoundStates.jl @@ -0,0 +1,59 @@ +#!/usr/bin/env julia + +"""Find eigenstates and eigenenergies of central potential problems""" + +using LinearAlgebra +using Plots + +N = 5000 # number of lattice points +L = 20.0 # r runs from 0 to L +dr = L / N + +D = zeros(N, N) # discrete radial 2nd derivative operator +V = zeros(N, N) # potential + +for i in 1:N + D[i, i] = -2.0 +end + +for i in 1:N-1 + D[i, i+1] = 1.0 + D[i+1, i] = 1.0 +end + +#println("\nLattice Laplacian operator") +#println(D) + +function potential(r, ℓ = 0) + """ The potential energy""" + #return 0.5 * ell * (ℓ+1.0) * pow(r, -2.0) # V=0: Free particle in spherical coordinates + + return -1.0/r + 0.5 * ℓ * (ℓ+1.0) * r^(-2.0) # Hydrogen atom + + #return -r^(-1.1) + 0.5 * ℓ * (ℓ+1.0) * r^(-2.0) # modified Coulomb potential +end + +for i in 1:N + r = (i + 0.5) * dr # radial coordinates of lattice points + V[i, i] = potential(r, 0) +end + +H = -0.5 * dr^(-2.0) * D + V # Hamiltonian. Here m = hbar = 1 + +#println("\nMatrix elements of Hamiltonian = ") +#println(H) + +e, v = eigen(H) # diagonalize Hamiltonian + +println("\nGround state energy = ", e[1]) +println("\n1st excited state energy = ", e[2]) +println("\n2nd excited state energy = ", e[3]) +println("\n3rd excited state energy = ", e[4]) +println("\n4th excited state energy = ", e[5]) + + +plot(potential) + + +plot(v[:,1]) +#plot(v[:,2]) diff --git a/hw7/mvsb.svg b/hw7/mvsb.svg new file mode 100644 index 0000000..9e4afbb --- /dev/null +++ b/hw7/mvsb.svg @@ -0,0 +1,55 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw7/mvsb2.svg b/hw7/mvsb2.svg new file mode 100644 index 0000000..9415495 --- /dev/null +++ b/hw7/mvsb2.svg @@ -0,0 +1,46 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw7/tightbinding (1).jl b/hw7/tightbinding (1).jl new file mode 100644 index 0000000..0d1cc8c --- /dev/null +++ b/hw7/tightbinding (1).jl @@ -0,0 +1,24 @@ +using LinearAlgebra + +# Define the tight binding Hamiltonian +function tight_binding_hamiltonian(n_sites::Int64, t::Float64) + H = zeros(n_sites, n_sites) + for i in 1:(n_sites-1) + H[i, i+1] = t + H[i+1, i] = t + end + H +end + +# Find the eigenenergies of the Hamiltonian +function eigenenergies(H::Array{Float64, 2}) + eigvals(H) +end + +n_sites = 10 # number of sites in the chain +t = 1.0 # hopping parameter +H = tight_binding_hamiltonian(n_sites, t) +energies = eigenenergies(H) +println(energies) + + -- cgit v1.2.3-70-g09d2