using Plots gr() ##ENV["PLOTS_TEST"] = "true" using SpecialFunctions using Distributions using SymPy lstirling(s) = s*log(s) - s - log(s)/2 + log(√(2π)) s = range(0.1, 10, length=400) plot(size=(500, 300), legend=:topleft) plot!(s, lgamma.(s), label="log Gamma(s)", lw=2) plot!(s, lstirling.(s), label="Stirling approximation", ls=:dash, lw=2) k, θ = 100, 1/80 gdist = Gamma(k, θ) ndist = Normal(k*θ, √k*θ) x = range(0, 2.7, length=400) plot(size=(500, 300)) title!("Normal approximation of Gamma distribution", titlefontsize=12) plot!(legend=:topright) plot!(x, pdf.(gdist, x), label="Gamma($k, $θ)", lw=2) plot!(x, pdf.(ndist, x), label="Normal approximation", ls=:dash, lw=2) a, b = 40, 60 bdist = Beta(a, b) ndist = Normal(a/(a+b), √(a*b/((a+b)^2*(a+b+1)))) x = range(0, 1, length=400) plot(size=(500, 300)) title!("Normal approximation of Beta distribution", titlefontsize=12) plot!(legend=:topright) plot!(x, pdf.(bdist, x), label="Beta($a, $b)", lw=2) plot!(x, pdf.(ndist, x), label="Normal approximation", ls=:dash, lw=2) k, θ, x = symbols("k θ x", positive=true) sol = simplify(integrate(exp(-x/θ)*x^(k-1), (x, 0, oo))) solstr = replace(sympy.latex(sol), "θ"=>"\\theta") display("text/html", raw"$$\int_0^\infty e^{-x/\theta} x^{k-1}\,dx = " * solstr * raw"$$"); f(s) = max(min(log(abs(zeta(s))), 3), -5) x = range(0, 1, length=100) y = range(10, 45, length=400) s = @. x + im*y' plot(size=(750, 180)) title!("log(abs(zeta(x+iy)))", titlefontsize=12) heatmap!(y, x, f.(s), color=:thermal) xlabel!("y") ylabel!("x") hline!([0.5], ls=:dot, color=:cyan, label="") f(s) = abs(zeta(s)) x = 1/2 y = range(10, 45, length=1000) s = @. x + im*y plot(size=(720, 200)) title!("abs(zeta(0.5+iy)))", titlefontsize=12) plot!(y, f.(s), label="") xlabel!("y") x = 1/2 y = range(0, 45, length=1000) s = @. x + im*y w = zeta.(s) plot(size=(500,400), legend=:topleft) title!("zeta(0.5+it) for 0 < t < 45", titlefontsize=12) plot!(real(w), imag(w), label="zeta(0.5+it)") scatter!([0], [0], markersize=2, label="") xlabel!("x") ylabel!("y") annotate!([(-1.25, 0.1, "zeta(0.5)", 9)]) annotate!([(2.55, 1.7, "zeta(0.5+45i)", 9)]) x = 0.6 y = range(0, 45, length=1000) s = @. x + im*y w = zeta.(s) plot(size=(500,400), legend=:topleft) title!("zeta(0.6+it) for 0 < t < 45", titlefontsize=12) plot!(real(w), imag(w), label="zeta(0.6+it)") scatter!([0], [0], markersize=2, label="") xlabel!("x") ylabel!("y") annotate!([(-1.6, 0.1, "zeta(0.6)", 9)]) annotate!([(2.55, 1.4, "zeta(0.6+45i)", 9)])