{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Homework 2 Solutions: Integrable Singularities and Monte Carlo Integration\n", "### Gabe Schumm\n", "### Due October 6th, 2021" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 1) Integrable Singularities" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Plots, Statistics, DelimitedFiles, LsqFit, LaTeXStrings, Printf\n", "#Plots for plotting\n", "#Statistics for mean\n", "#DelimitedFiles to write data (not used here)\n", "#LsqFit to fit to do straight line fit\n", "#LaTeXStrings for nice plot font\n", "#Printf for string formating" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I will first define a module that will contain the Integrable Singularities structure `IntSing`. The reason I do this is that if we want to edit our struct in any way (add a new field for example) we can do so in the structure defined within the module with not issues. If the structure was defined globally, it will be uneditable and you would have to restart your Kernel in order to adjust it.\n", "\n", "I add `mutable` in front of `struct` so that the values of the fields for a certain instance can be changed after its initially defined. For this function. This isn't necessary but in other cases, this will be helpful." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "module IntSings\n", "\n", "mutable struct IntSing #define the structure IntSing that will contain the specifications of the integration\n", " ϵ::Float64\n", " α::Float64\n", " \n", " x0::Float64\n", " xf::Float64\n", "\n", " N0::Int64\n", " N::Int64\n", " n::Int64\n", "\n", " h::Float64\n", " \n", " ifxn::Any #function we will use for integration (either first or second order)\n", "end\n", "\n", "\n", "#define first order integration function within module\n", "function int1(self)\n", " function f(x) return 1/(self.ϵ + x)^self.α end\n", " \n", " I = 1.5*(f(self.x0+self.h) + f(self.xf-self.h)) # add end points\n", " \n", " for i=2:self.N-2 # n must be even\n", " I += f(self.x0+ (self.h* i))\n", " end\n", " return I * self.h\n", "end\n", " \n", "#define second order integration function within module\n", "function int2(self)\n", " function f(x) return 1/(self.ϵ + x)^self.α end\n", " \n", " I = ((23/12)*(f(self.x0+self.h) + f(self.xf-self.h)) +\n", " (7/12)*(f(self.x0+(2*self.h)) + f(self.xf-(2*self.h))))# add end points\n", " \n", " for i=3:self.N-3 # n must be even\n", " I += f(self.x0+ (self.h* i))\n", " end\n", " return I * self.h\n", "end\n", "\n", "\n", "#function with same name as structure, used to initialize structure instance and set constants\n", "function IntSing(ϵ, α, n, ifxn_id)\n", " \n", " #these are the same for all runs\n", " x0 = 0\n", " xf = 1\n", " N0 = 10\n", " \n", " N = N0 * 2^n\n", " \n", " #if user input ifxn_id=1, use first order, if 2 use second order\n", " if ifxn_id == 1\n", " ifxn = int1\n", " elseif ifxn_id == 2\n", " ifxn = int2\n", " else\n", " error(\"Invalid ifxn_id.\")\n", " end\n", "\n", " \n", " h = (xf - x0)/N\n", " \n", " #return instance of IntSing with user specified and default variables\n", " return IntSing(ϵ, α,\n", " x0, xf,\n", " N0, N, n,\n", " h,\n", " ifxn)\n", " \n", "\n", "end\n", "\n", "end\n", "\n", "#define IntSing as IntSing from module IntSings in global scope\n", "IntSing = IntSings.IntSing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part a" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "linear (generic function with 1 method)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "linear(x, p) = p[1] .+ x .* p[2] #function used in curve fit (y=mx+b)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hs = []\n", "I1 = [] #first order\n", "I2 = [] #second order\n", "\n", "ϵ = 1\n", "α = 0.5\n", "for n=0:11\n", " intsing1 = IntSing(ϵ, α, n, 1) #define instance with first order integral funciton\n", " push!(I1, intsing1.ifxn(intsing1))\n", " \n", " intsing2 = IntSing(ϵ, α, n, 2) #define instance with second order integral funciton\n", " push!(I2, intsing2.ifxn(intsing2))\n", " \n", " push!(hs, intsing2.h)\n", "end\n", "\n", "I_exact = 2^(1.5) - 2 #exact formula for integral\n", "\n", "#fit straight line to log-log data with arbitrary initial parameters (.5, .5)\n", "fit_1 = curve_fit(linear, log.(hs), log.(abs.(I1 .- I_exact)), [.5,.5]) \n", "fit_2 = curve_fit(linear, log.(hs), log.(abs.(I2 .- I_exact)), [.5,.5]);\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "scatter(log.(hs), log.(abs.(I1 .- I_exact)), label =\"First Order\", c=\"blue\", legend=:right,\n", " xlabel = L\"\\log(h)\", ylabel =L\"\\Delta_k(N)\")\n", "plot!(log.(hs), linear((log.(hs)), fit_1.param), label =\"\", c=\"blue\")\n", "\n", "\n", "string_1 = latexstring( L\"\\gamma_1=\",@sprintf(\"%.3f\", fit_1.param[2]))\n", "\n", "\n", "scatter!(log.(hs), log.(abs.(I2 .- I_exact)), label =\"Second Order\", c=\"red\")\n", "plot!(log.(hs), linear((log.(hs)), fit_2.param), label =\"\", c=\"red\")\n", "\n", "string_2 = latexstring( L\"\\gamma_2=\",@sprintf(\"%.3f\", fit_2.param[2]))\n", "annotate!([(-7, -10, string_1), (-5, -25, string_2)])\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "The exponents agree with the expected values." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part b" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "ϵ = 0\n", "α = .5\n", "\n", "hs = []\n", "I1 = []\n", "I2 = []\n", "for n=0:20\n", " #println(n)\n", " intsing1 = IntSing(ϵ, α, n, 1)\n", " push!(I1, intsing1.ifxn(intsing1))\n", " \n", " intsing2 = IntSing(ϵ, α, n, 2)\n", " push!(I2, intsing2.ifxn(intsing2))\n", " \n", " push!(hs, intsing2.h)\n", "end\n", "\n", "I_exact = 1/(-α + 1)\n", "\n", "\n", "fit_1 = curve_fit(linear, log.(hs), log.(abs.(I1 .- I_exact)), [.5,.5])\n", "fit_2 = curve_fit(linear, log.(hs), log.(abs.(I2 .- I_exact)), [.5,.5]);\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "scatter(log.(hs), log.(abs.(I1 .- I_exact)), label =\"First Order\", c=\"blue\", legend=:right,\n", " xlabel = L\"\\log(h)\", ylabel =L\"\\Delta_k(N)\")\n", "plot!(log.(hs), linear(log.(hs), fit_1.param), label =\"\", c=\"blue\")\n", "\n", "\n", "string_1 = latexstring( L\"\\gamma_1=\",@sprintf(\"%.5f\", fit_1.param[2]))\n", "\n", "\n", "scatter!(log.(hs), log.(abs.(I2 .- I_exact)), label =\"Second Order\", c=\"red\")\n", "plot!(log.(hs), linear(log.(hs), fit_2.param), label =\"\", c=\"red\", title = L\"\\alpha = 0.5\")\n", "\n", "string_2 = latexstring( L\"\\gamma_2=\",@sprintf(\"%.5f\", fit_2.param[2]))\n", "\n", "annotate!([(-10, -3, string_1), (-7, -7, string_2)])\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ϵ = 0\n", "α = .75\n", "\n", "hs = []\n", "I1 = []\n", "I2 = []\n", "for n=0:20\n", " #println(n)\n", " intsing1 = IntSing(ϵ, α, n, 1)\n", " push!(I1, intsing1.ifxn(intsing1))\n", " \n", " intsing2 = IntSing(ϵ, α, n, 2)\n", " push!(I2, intsing2.ifxn(intsing2))\n", " \n", " push!(hs, intsing2.h)\n", "end\n", "\n", "I_exact = 1/(-α + 1)\n", "\n", "\n", "fit_1 = curve_fit(linear, log.(hs), log.(abs.(I1 .- I_exact)), [.5,.5])\n", "fit_2 = curve_fit(linear, log.(hs), log.(abs.(I2 .- I_exact)), [.5,.5]);\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "scatter(log.(hs), log.(abs.(I1 .- I_exact)), label =\"First Order\", c=\"blue\", legend=:right,\n", " xlabel = L\"\\log(h)\", ylabel =L\"\\Delta_k(N)\")\n", "plot!(log.(hs), linear(log.(hs), fit_1.param), label =\"\", c=\"blue\")\n", "\n", "\n", "string_1 = latexstring( L\"\\gamma_1=\",@sprintf(\"%.5f\", fit_1.param[2]))\n", "\n", "\n", "scatter!(log.(hs), log.(abs.(I2 .- I_exact)), label =\"Second Order\", c=\"red\")\n", "plot!(log.(hs), linear(log.(hs), fit_2.param), label =\"\", c=\"red\", title = L\"\\alpha = 0.75\")\n", "\n", "string_2 = latexstring( L\"\\gamma_2=\",@sprintf(\"%.5f\", fit_2.param[2]))\n", "\n", "annotate!([(-12, -0.5, string_1), (-7, -2, string_2)])\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In both cases, the exponent $\\gamma$ is\n", "$$ \\gamma = 1-\\alpha$$\n", "\n", "For the open interval formula, the error is dominated by the singular end point at $x=0$:\n", "\n", "$$\n", "\\Delta I = |I_k(N) - I| \\approx \\int_0^h f(x) dx\n", "$$\n", "where the integral on the right-hand side is the *exlcuded* contribution to the integral in our forumals. We can solve this explictly to calculate the $h$ dependance.\n", "$$\n", "\\int_0^h f(x) dx = \\int_0^h \\frac{1}{x^\\alpha} dx = \\frac{x^{1-\\alpha}}{1-\\alpha}\\Big|_0^h = \\frac{h^{1-\\alpha}}{1-\\alpha}\\propto h^{1-\\alpha}\n", "$$\n", "\n", "Thus, the dominant contribution to the error from excluding the singular point scales as $h^{1-\\alpha}$, confirming the results of the above plots. Note that there is still a $\\mathcal{O}\\left(h^2\\right) or \\mathcal{O}\\left(h^3\\right)$ error contribution, but the $\\mathcal{O}\\left(h^{1-\\alpha}\\right)$ contribution dominates." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "α = .5\n", "ϵ = 1e-5\n", "hs = []\n", "I1 = []\n", "I2 = []\n", "for n=0:21\n", " #println(n)\n", " intsing1 = IntSing(ϵ, α,n, 1)\n", " push!(I1, intsing1.ifxn(intsing1))\n", " \n", " intsing2 = IntSing(ϵ, α, n, 2)\n", " push!(I2, intsing2.ifxn(intsing2))\n", " \n", " push!(hs, intsing2.h)\n", "end\n", "\n", "I_exact = ((1+ϵ)^(1-α)) / (1-α) - ((ϵ)^(1-α)) / (1-α)\n", "\n", "\n", "fit_1a = curve_fit(linear, log.(hs)[1:4], log.(abs.(I1 .- I_exact))[1:4], [.5,.5])\n", "fit_2a = curve_fit(linear, log.(hs)[1:4], log.(abs.(I2 .- I_exact))[1:4], [.5,.5]);\n", "\n", "\n", "fit_1b = curve_fit(linear, log.(hs)[18:end], log.(abs.(I1 .- I_exact))[18:end], [.5,.5])\n", "fit_2b = curve_fit(linear, log.(hs)[18:end], log.(abs.(I2 .- I_exact))[18:end], [.5,.5]);\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "scatter(log.(hs), log.(abs.(I1 .- I_exact)), label =\"First Order\", c=\"blue\", legend=:bottomright,\n", " xlabel = L\"\\log(h)\", ylabel =L\"\\Delta_k(N)\")\n", "plot!(log.(hs)[1:13], linear(log.(hs)[1:13], fit_1a.param), label =\"\", c=\"blue\")\n", "plot!(log.(hs)[12:end], linear(log.(hs)[12:end], fit_1b.param), label =\"\", c=\"blue\")\n", "\n", "\n", "string_1a = latexstring( L\"\\gamma_1=\",@sprintf(\"%.5f\", fit_1a.param[2]))\n", "string_1b = latexstring( L\"\\gamma_1=\",@sprintf(\"%.5f\", fit_1b.param[2]))\n", "\n", "\n", "scatter!(log.(hs), log.(abs.(I2 .- I_exact)), label =\"Second Order\", c=\"red\")\n", "plot!(log.(hs)[1:13], linear(log.(hs)[1:13], fit_2a.param), label =\"\", c=\"red\")\n", "plot!(log.(hs)[12:end], linear(log.(hs)[12:end], fit_2b.param), label =\"\", c=\"red\")\n", "\n", "\n", "\n", "string_2a = latexstring( L\"\\gamma_2=\",@sprintf(\"%.5f\", fit_2a.param[2]))\n", "string_2b = latexstring( L\"\\gamma_2=\",@sprintf(\"%.5f\", fit_2b.param[2]))\n", "\n", "\n", "annotate!([(-6, -8, string_1a), (-6, -10, string_2a),\n", " (-14, -4, string_1b), (-14, -6, string_2b)])\n", "\n", "plot!([log(ϵ), log(ϵ)], [0, -23], linestyle=:dash, label = L\"x=\\log(\\epsilon)\")\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We do indeed see a \"cross-over\" region, where the scaling of the error switches from $\\mathcal{O}\\left(h^{1-\\alpha}\\right)$ to $\\mathcal{O}\\left(h^2\\right) or \\mathcal{O}\\left(h^3\\right)$. \n", "\n", "When $h$ is larger than $\\epsilon$ (to the right of the dotted black line), we see the $\\mathcal{O}\\left(h^{1-\\alpha}\\right)$ scaling. In this regime, the step size is larger than the deviation from the singular form, so it's as if we don't have the resolution to determine that $f(x)$ is not divergent at $x=0$. \n", "\n", "When $h$ is less than $\\epsilon$ (to the left of the dotted black line), we see the $\\mathcal{O}\\left(h^2\\right) and \\mathcal{O}\\left(h^3\\right)$ scaling. In this regime, the step size is small than the deviation from the singular form, so we can resolve the non-divergence of $f(x)$ is at $x=0$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 1) Monte Carlo Integration\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "using CSV, DataFrames" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulation Functions" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "analyze_data (generic function with 1 method)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Sphere structure containing\n", "r_1: radius of sphere\n", "r_2: radius of cylinder\n", "rho_1: density of sphere\n", "rho_2: density of cylinder\n", "z_cyl: positive z position of top of cylinder (h/2)\n", "L: length of box enclosing sphere (2* r_1)\n", "bin_outfile: path for bin average results\n", "res_outfile: path for final average and error bar results\n", "\"\"\"\n", "struct Sphere\n", " r_1::Float64\n", " r_2::Float64\n", "\n", " rho_1::Float64\n", " rho_2::Float64\n", " \n", " z_cyl::Float64\n", " \n", " L::Float64\n", "\n", " bin_outfile::String\n", " res_outfile::String\n", "end\n", "\n", "\n", "function Sphere(r_1, r_2, rho_1, rho_2, bin_outfile, res_outfile)\n", " \n", " z_cyl = sqrt(r_1^2 - r_2^2) #z position of clyinder found using Pythagorean theorem\n", " \n", " sphere = Sphere(r_1, r_2, rho_1, rho_2, z_cyl, 2*r_1, bin_outfile, res_outfile)\n", " \n", "end\n", "\n", "#Function that returns whether or not point (x, y, z) is in the cylinder\n", "function in_cylinder(sphere, x, y, z)\n", " ρ = (x^2 + y^2) ^ (1/2) #radius in cylindrical coordiantes\n", " \n", " if ρ <= sphere.r_2 && abs(z) <= sphere.z_cyl #condition that point (x, y, z) is in the cylinder\n", " return true\n", " else\n", " return false\n", " end\n", "end\n", "\n", "#value of moment of inertia integrand at point x, y, z, for moment of interatia along given axis\n", "function integrand(sphere, x, y, z, axis)\n", " \n", " #perpendicular distance from axis\n", " if axis == \"x\"\n", " r_perp2 = y^2 + z^2\n", " elseif axis == \"z\"\n", " r_perp2 = x^2 + y^2\n", " else\n", " throw(\"Invalid axis of rotation\")\n", " end\n", " \n", " \n", " if in_cylinder(sphere, x, y, z)\n", " return r_perp2 * sphere.rho_2\n", " else\n", " return r_perp2 * sphere.rho_1\n", " end\n", "end\n", "\n", "#Integrate over sphere object using npt points per simulation for nbin simulations \n", "function integrate(sphere, npt, nbin)\n", " Ix = []\n", " Iz = []\n", " \n", " \n", " #open empty bin average output file\n", " #overwrite if file already exists\n", " open(sphere.bin_outfile, \"w\") do f\n", "\n", " for nb = 1:nbin\n", " \n", " #generate all npt (x, y, z) samples for the simulation\n", " #(x, y, z) = 3 random numbers between -L and L\n", " samples = rand(npt, 3) .* (2*sphere.r_1) .- (sphere.r_1)\n", " \n", " #print status of simulation\n", " if nb % 1000 == 0\n", " println(\"Bin number: \",nb)\n", " end\n", " \n", " #initialize monent of intertias for each simulations\n", " Ix_bin = 0 \n", " Iz_bin = 0 \n", " for np = 1:npt\n", " #select np'th set of (x, y, z) values\n", " x, y, z = samples[np, :]\n", "\n", " # if (x, y, z) is not in sphere, return zero for value of integrand\n", " if (x^2 + y^2 + z^2) > sphere.r_1^2\n", " Ix_bin += 0\n", " Iz_bin += 0\n", " #else add 1/npt * value of integrand at (x, y, z), for each monent of inertia\n", " #(dividing by 1/npt in this step takes the avearge over all samples)\n", " else\n", " Ix_bin += integrand(sphere, x, y, z, \"x\") / npt\n", " Iz_bin += integrand(sphere, x, y, z, \"z\") / npt\n", " end\n", " end\n", " \n", " #write simulation number, I_z, and I_x to file\n", " #multiply integrand average by box volume sphere.L^3 to get moment of inertia\n", " println(f, nb, \" \", Iz_bin*(sphere.L^3), \" \", Ix_bin*(sphere.L^3))\n", "\n", " end\n", " end\n", " \n", "end\n", "\n", "function V_sphere(sphere, npt, nbin)\n", " V = []\n", " \n", " for nb in 1:nbin\n", " V_bin = 0 \n", " for np in 1:npt\n", " x, y, z = rand(3) .* (sphere.L) .- (sphere.L/2)\n", " \n", " if (x^2 + y^2 + z^2) > sphere.r_1\n", " V_bin += 0\n", " else\n", " V_bin += 1/npt\n", " end\n", " end\n", " push!(V, V_bin)\n", " end\n", " return V\n", "end\n", "\n", "function analyze_data(sphere)\n", " #load bin averages as DataFrame\n", " df = CSV.read(sphere.bin_outfile, delim = \" \", header=false, DataFrame)\n", " \n", " \n", " nbins = [50, 500, 5000]\n", " \n", " #calclate statistics for each bin number\n", " #write statistics to result outfile\n", " open(sphere.res_outfile, \"w\") do f\n", " for nbin in nbins\n", " Iz_sample = df[1:nbin, \"Column2\"]\n", " Ix_sample = df[1:nbin, \"Column3\"]\n", "\n", " σ_z = sqrt(1/(nbin*(nbin-1)) * sum(Iz_sample .^2 .- mean(Iz_sample)^2))\n", " σ_x = sqrt(1/(nbin*(nbin-1)) * sum(Ix_sample .^2 .- mean(Ix_sample)^2))\n", " \n", " println(f, nbin, \" \", mean(Iz_sample), \" \", σ_z, \" \", mean(Ix_sample), \" \", σ_x) \n", "\n", " end\n", " end\n", " println(\"Statistics exported\")\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initializing Sphere Object" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "s = Sphere(.05, .01, 8930, 19320, \"bin.dat\", \"res.dat\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running 5,000 simulations" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bin number: 1000\n", "Bin number: 2000\n", "Bin number: 3000\n", "Bin number: 4000\n", "Bin number: 5000\n", "334.361077 seconds (5.00 G allocations: 745.070 GiB, 6.49% gc time, 0.03% compilation time)\n" ] } ], "source": [ "@time integrate(s, 10^6, 5000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculating Statistics for $n_{bin}$ = 50, 500, 5,000" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Statistics exported\n" ] } ], "source": [ "analyze_data(s)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Comparing Results to Exact Value of $I_z, I_x$" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/html": [ "

3 rows × 5 columns

Column1Column2Column3Column4Column5
Int64Float64Float64Float64Float64
1500.004692728.14559e-70.004939331.04907e-6
25000.004691922.61004e-70.00493963.18356e-7
350000.004691788.74882e-80.004939659.62373e-8
" ], "text/latex": [ "\\begin{tabular}{r|ccccc}\n", "\t& Column1 & Column2 & Column3 & Column4 & Column5\\\\\n", "\t\\hline\n", "\t& Int64 & Float64 & Float64 & Float64 & Float64\\\\\n", "\t\\hline\n", "\t1 & 50 & 0.00469272 & 8.14559e-7 & 0.00493933 & 1.04907e-6 \\\\\n", "\t2 & 500 & 0.00469192 & 2.61004e-7 & 0.0049396 & 3.18356e-7 \\\\\n", "\t3 & 5000 & 0.00469178 & 8.74882e-8 & 0.00493965 & 9.62373e-8 \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "\u001b[1m3×5 DataFrame\u001b[0m\n", "\u001b[1m Row \u001b[0m│\u001b[1m Column1 \u001b[0m\u001b[1m Column2 \u001b[0m\u001b[1m Column3 \u001b[0m\u001b[1m Column4 \u001b[0m\u001b[1m Column5 \u001b[0m\n", "\u001b[1m \u001b[0m│\u001b[90m Int64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\n", "─────┼─────────────────────────────────────────────────────────\n", " 1 │ 50 0.00469272 8.14559e-7 0.00493933 1.04907e-6\n", " 2 │ 500 0.00469192 2.61004e-7 0.0049396 3.18356e-7\n", " 3 │ 5000 0.00469178 8.74882e-8 0.00493965 9.62373e-8" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = CSV.read(\"res.dat\", delim = \" \", DataFrame, header=false)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact Error: 4.913184445138469e-8\n", "Monte Carlo σ_z: 8.74882243583041e-8\n" ] } ], "source": [ "#Exact I_z\n", "I_cyl_rho2z = 19320 * pi *.01^4 * sqrt(.05^2 - .01^2)\n", "I_cyl_rho1z = 8930 * pi *.01^4 * sqrt(.05^2 - .01^2)\n", "\n", "I_sph_rho1 = 8930 * (8 *pi /15) *.05^5 \n", " \n", "error_z = abs((I_sph_rho1 - I_cyl_rho1z + I_cyl_rho2z) - df[3, \"Column2\"])\n", "\n", "println(\"Exact Error: \",error_z)\n", "println(\"Monte Carlo σ_z: \",df[3, \"Column3\"])" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact Error: -5.983015735535474e-8\n", "Monte Carlo σ_x: 9.623728119689116e-8\n" ] } ], "source": [ "#Exact I_x\n", "I_cyl_rho2x = 19320 * (pi/6) *.01^2 * sqrt(.05^2 - .01^2) * ((4*.05^2) - .01^2)\n", "I_cyl_rho1x = 8930 * (pi/6) *.01^2 * sqrt(.05^2 - .01^2) * ((4*.05^2) - .01^2)\n", "\n", "\n", "error_x = (I_sph_rho1 - I_cyl_rho1x + I_cyl_rho2x) - df[3, \"Column4\"]\n", "\n", "println(\"Exact Error: \",error_x)\n", "println(\"Monte Carlo σ_x: \",df[3, \"Column5\"])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.6.2", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.2" } }, "nbformat": 4, "nbformat_minor": 2 }