143 lines
3.3 KiB
Julia
143 lines
3.3 KiB
Julia
|
using AbstractPlotting
|
|||
|
|
using AbstractPlotting.MakieLayout
|
|||
|
|
using GLMakie
|
|||
|
|
|
|||
|
|
using QuadGK
|
|||
|
|
|
|||
|
|
function riemann(f::Function, a::Real, b::Real, n::Int; method="right")
|
|||
|
|
if method == "right"
|
|||
|
|
meth = f -> (lr -> begin l,r = lr; f(r) * (r-l) end)
|
|||
|
|
elseif method == "left"
|
|||
|
|
meth = f -> (lr -> begin l,r = lr; f(l) * (r-l) end)
|
|||
|
|
elseif method == "trapezoid"
|
|||
|
|
meth = f -> (lr -> begin l,r = lr; (1/2) * (f(l) + f(r)) * (r-l) end)
|
|||
|
|
elseif method == "simpsons"
|
|||
|
|
meth = f -> (lr -> begin l,r=lr; (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l) end)
|
|||
|
|
end
|
|||
|
|
|
|||
|
|
xs = a .+ (0:n) * (b-a)/n
|
|||
|
|
|
|||
|
|
sum(meth(f), zip(xs[1:end-1], xs[2:end]))
|
|||
|
|
end
|
|||
|
|
|
|||
|
|
|
|||
|
|
"""
|
|||
|
|
integration(f)
|
|||
|
|
|
|||
|
|
Show graphically the left Riemann approximation, the trapezoid approximation, and Simpson's approximation to the integral of `f` over [-1,1].
|
|||
|
|
|
|||
|
|
Assumes `f` is non-negative.
|
|||
|
|
"""
|
|||
|
|
function integration(f=nothing)
|
|||
|
|
if f == nothing
|
|||
|
|
f = x -> x^2*exp(-x/3)
|
|||
|
|
end
|
|||
|
|
|
|||
|
|
a, b = -1, 1
|
|||
|
|
|
|||
|
|
|
|||
|
|
function left_riemann_pts(n)
|
|||
|
|
xs = range(a, b, length=n+1)
|
|||
|
|
pts = Point2f0[(xs[1], 0)]
|
|||
|
|
for i in 1:n
|
|||
|
|
xᵢ, xᵢ₊₁ = xs[i], xs[i+1]
|
|||
|
|
fᵢ = f(xᵢ)
|
|||
|
|
push!(pts, (xᵢ, fᵢ))
|
|||
|
|
push!(pts, (xᵢ₊₁, fᵢ))
|
|||
|
|
end
|
|||
|
|
push!(pts, (xs[end], 0))
|
|||
|
|
pts
|
|||
|
|
end
|
|||
|
|
|
|||
|
|
|
|||
|
|
function trapezoid_pts(n)
|
|||
|
|
xs = range(a, b, length=n+1)
|
|||
|
|
pts = Point2f0[(xs[1], 0)]
|
|||
|
|
for i in 1:n
|
|||
|
|
xᵢ, xᵢ₊₁ = xs[i], xs[i+1]
|
|||
|
|
fᵢ = f(xᵢ)
|
|||
|
|
push!(pts, (xᵢ, f(xᵢ)))
|
|||
|
|
end
|
|||
|
|
push!(pts, (xs[end], f(xs[end])))
|
|||
|
|
push!(pts, (xs[end], 0))
|
|||
|
|
pts
|
|||
|
|
end
|
|||
|
|
|
|||
|
|
function simpsons_pts(n)
|
|||
|
|
xs = range(a, b, length=n+1)
|
|||
|
|
pts = Point2f0[(xs[1], 0), (xs[1], f(xs[1]))]
|
|||
|
|
|
|||
|
|
for i in 1:n
|
|||
|
|
xi, xi1 = xs[i], xs[i+1]
|
|||
|
|
m = xi/2 + xi1/2
|
|||
|
|
p = x -> f(xi)*(x-m)*(x - xi1)/(xi-m)/(xi-xi1) + f(m) * (x-xi)*(x-xi1)/(m-xi)/(m-xi1) + f(xi1) * (x-xi)*(x-m) / (xi1-xi) / (xi1-m)
|
|||
|
|
|
|||
|
|
xs′ = range(xi, xi1, length=10)
|
|||
|
|
for j in 2:10
|
|||
|
|
x = xs′[j]
|
|||
|
|
push!(pts, (x, p(x)))
|
|||
|
|
end
|
|||
|
|
end
|
|||
|
|
push!(pts, (xs[end], 0))
|
|||
|
|
pts
|
|||
|
|
end
|
|||
|
|
|
|||
|
|
|
|||
|
|
# where we lay our scene:
|
|||
|
|
scene, layout = layoutscene()
|
|||
|
|
layout.halign = :left
|
|||
|
|
layout.valign = :top
|
|||
|
|
|
|||
|
|
|
|||
|
|
p1 = layout[1,1] = LScene(scene)
|
|||
|
|
p2 = layout[1,2] = LScene(scene)
|
|||
|
|
p3 = layout[1,3] = LScene(scene)
|
|||
|
|
|
|||
|
|
n = layout[2,1] = LSlider(scene, range=2:25, startvalue=2)
|
|||
|
|
output = layout[2, 2:3] = LText(scene, "...")
|
|||
|
|
|
|||
|
|
lpts = lift(n.value) do n
|
|||
|
|
left_riemann_pts(n)
|
|||
|
|
end
|
|||
|
|
|
|||
|
|
poly!(p1.scene, lpts, color=:gray75)
|
|||
|
|
lines!(p1.scene, a..b, f, color=:black, strokewidth=10)
|
|||
|
|
title(p1.scene, "Left Riemann")
|
|||
|
|
|
|||
|
|
|
|||
|
|
tpts = lift(n.value) do n
|
|||
|
|
trapezoid_pts(n)
|
|||
|
|
end
|
|||
|
|
poly!(p2.scene, tpts, color=:gray75)
|
|||
|
|
lines!(p2.scene, a..b, f, color=:black, strokewidth=10)
|
|||
|
|
title(p2.scene, "Trapezoid")
|
|||
|
|
|
|||
|
|
spts = lift(n.value) do n
|
|||
|
|
simpsons_pts(n)
|
|||
|
|
end
|
|||
|
|
poly!(p3.scene, spts, color=:gray75)
|
|||
|
|
lines!(p3.scene, a..b, f, color=:black, strokewidth=10)
|
|||
|
|
title(p3.scene, "Simpson's")
|
|||
|
|
|
|||
|
|
on(n.value) do n
|
|||
|
|
actual,err = quadgk(f, -1, 1)
|
|||
|
|
lrr = riemann(f, -1, 1, n, method="left")
|
|||
|
|
trap = riemann(f, -1, 1, n, method="trapezoid")
|
|||
|
|
simp = riemann(f, -1, 1, n, method="simpsons")
|
|||
|
|
|
|||
|
|
Δleft = round(abs(lrr - actual), digits=8)
|
|||
|
|
Δtrap = round(abs(trap - actual), digits=8)
|
|||
|
|
Δsimp = round(abs(simp - actual), digits=8)
|
|||
|
|
|
|||
|
|
txt = "Riemann: $Δleft, Trapezoid: $Δtrap, Simpson's: $Δsimp"
|
|||
|
|
output.text[] = txt
|
|||
|
|
end
|
|||
|
|
n.value[] = n.value[]
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
scene
|
|||
|
|
end
|