diff --git a/run.jl b/run.jl index 508df9c..a885f86 100644 --- a/run.jl +++ b/run.jl @@ -1,4 +1,6 @@ -include("src/HW_int.jl") # load our code -HWintegration.runall() \ No newline at end of file +include("src/HWintegration.jl") # load our code +HWintegration.runall() + +plot(p...) diff --git a/src/HWintegration.jl b/src/HWintegration.jl index 68ef546..1ad6fb3 100644 --- a/src/HWintegration.jl +++ b/src/HWintegration.jl @@ -1,7 +1,7 @@ module HWintegration - const A_SOL = -18 # enter your analytic solution. -18 is wrong. + const A_SOL = 4 # enter your analytic solution. # imports using FastGaussQuadrature @@ -15,51 +15,84 @@ module HWintegration # set random seed srand(12345) - + p=Any[] # demand function + q(x)=2/(sqrt(x)) + a=3/2 + b=5/2 + g(x)=a*q(a*x+b) - # gauss-legendre adjustment factors for map change + # Question 1 + # gauss-legendre adjustment factors + a=3/2 + b=5/2 # eqm condition for question 2 # makes a plot for questions 1 function plot_q1() - + plot(q,0,4,xlim=(0,4), ylim=(0,4),label="Demand") + hline!([q(1)],label="p=1") + hline!([q(4)],label="p=2") end function question_1b(n) - + println(sum(gausslegendre(n)[2][i]*g(gausslegendre(n)[1][i]) for i in 1:n )), + push!(p,plot(g,-1,1,xlim=(-1,1), ylim=(0,4),title="GaussLegendre",legend=false)) + z=[] + for i in 1:n + push!(z,g(gausslegendre(n)[1][i])) + end + scatter!(gausslegendre(n)[1],z,legend=false) end function question_1c(n) + MonteCarlo=rand(Uniform(1,4),n) + println(3*sum((1/n)*q(MonteCarlo[i]) for i in 1:n )), + push!(p,plot(q,0,4,xlim=(0,4), ylim=(0,4),title="MonteCarlo",legend=false)) + z=[] + for i in 1:n + push!(z,q(MonteCarlo[i])) + end + scatter!(MonteCarlo,z,legend=false) - end function question_1d(n) + QMC=ScaledSobolSeq(1,[1.0],[4.0]) + QuasiMonteCarlo=([next(QMC) for i in 1:n]) + QMCX=[QuasiMonteCarlo[i][1] for i in 1:n] + println(3*sum((1/n)*q(QMCX[i]) for i in 1:n )), + push!(p,plot(q,0,4,xlim=(0,4), ylim=(0,4),title="QuasiMonteCarlo",legend=false)) + z=[] + for i in 1:n + push!(z,q(QuasiMonteCarlo[i][1])) + end + scatter!(QMCX,z,legend=false) - end # question 2 function question_2a(n) + println("I could not answer this question") - end function question_2b(n) + println("I could not answer this question") + - end function question_2bonus(n) + println("I could not answer this question") + - end @@ -68,17 +101,18 @@ module HWintegration function runall() info("Running all of HWintegration") info("question 1:") - plot_q1() + plot_q1 for n in (10,15,20) info("============================") info("now showing results for n=$n") - # info("question 1b:") - # question_1b(n) # make sure your function prints some kind of result! - # info("question 1c:") - # question_1c(n) - # info("question 1d:") - # question_1d(n) - # println("") + info("question 1b:") + question_1b(n) # make sure your function prints some kind of result! + info("question 1c:") + question_1c(n) + info("question 1d:") + question_1d(n) + + println("") info("question 2a:") q2 = question_2a(n) println(q2) @@ -88,7 +122,9 @@ module HWintegration info("bonus question: Quasi monte carlo:") q2bo = question_2bonus(n) println(q2bo) + info("Creating Plot") end + plot(p...) end info("end of HWintegration")