# For CalculusWithJulia, is there Riemann' Midpoint method?

https://docs.juliahub.com/CalculusWithJulia/AZHbv/0.0.5/integrals/area.html

and I want to calculate Riemann sum for several functions, comparing the

1. Left
2. Right
3. Midpoint

methods. Is the midpoint method already available? If not then how to add it to the current Riemann Sum function?

1 Like

No, that function is for illustrative purposes only and doesnβt allow modification of the provided methods. You could certainly suggest some modification to the package, but for now, you might follow this simple enough pattern:

function Riemann(f, a, b, n; meth=(a,b) -> (a+b)/2)
xs = range(a, b, n+1)
pairs = zip(xs[begin:end-1], xs[begin+1:end]) # (xβ,xβ), β¦, (xβββ,xβ)
sum(f(meth(xs...)) for xs in pairs)
end


Pass in meth=(a,b) -> a for left sums, meth=(a,b)->b for rightβ¦

BTW, that link is out of date. Try https://calculuswithjulia.github.io for a more current version.

1 Like

Hi @j_verzani ,

I was using this function code:

function riemann(f, a, b, n; method="left")
xs = a:(b-a)/n:b
deltas = diff(xs)      # forms x2-x1, x3-x2, ..., xn-xn-1
if method == "left"
cs = xs[1:end-1]
elseif method == "right"
cs = xs[2:end]
elseif method == "midpoint"
cs = [(xs[i] + xs[i+1])/2 for i in 1:length(deltas)]
end
sum(f(cs[i]) * deltas[i] for i in 1:length(deltas))
end


and it can calculate the midpoint method correctly.

I think it is another version of your suggested function:

function Riemann(f, a, b, n; meth=(a,b) -> (a+b)/2)
xs = range(a, b, n+1)
pairs = zip(xs[begin:end-1], xs[begin+1:end]) # (xβ,xβ), β¦, (xβββ,xβ)
sum(f(meth(xs...)) for xs in pairs)
end


How to check which function works better? with btime?

If all the materials in (Calculus with Julia) printed to become a book I will look forward to it, it is really helping me to study better than in classes.

I have been looking at this function from the old blog of yours with midpoint method:

function riemann(f, a, b, n; method="left")
xs = a:(b-a)/n:b
deltas = diff(xs)      # forms x2-x1, x3-x2, ..., xn-xn-1
if method == "left"
cs = xs[1:end-1]
elseif method == "right"
cs = xs[2:end]
elseif method == "midpoint"
cs = [(xs[i] + xs[i+1])/2 for i in 1:length(deltas)]
end
sum(f(cs[i]) * deltas[i] for i in 1:length(deltas))
end


I wonder how to convert the midpoint method into the new function in the updated blog, so all methods will be available:

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 = range(a, b, n+1)
pairs = zip(xs[begin:end-1], xs[begin+1:end]) # (xβ,xβ), β¦, (xβββ,xβ)
sum(meth(f), pairs)

end


My inexperience resulting in this to make all methods available:

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)
xs = range(a, b, n+1)
pairs = zip(xs[begin:end-1], xs[begin+1:end]) # (xβ,xβ), β¦, (xβββ,xβ)
sum(meth(f), pairs)
elseif method == "left"
meth = f -> (lr -> begin l,r = lr; f(l) * (r-l) end)
xs = range(a, b, n+1)
pairs = zip(xs[begin:end-1], xs[begin+1:end]) # (xβ,xβ), β¦, (xβββ,xβ)
sum(meth(f), pairs)
elseif method == "trapezoid"
meth = f -> (lr -> begin l,r = lr; (1/2) * (f(l) + f(r)) * (r-l) end)
xs = range(a, b, n+1)
pairs = zip(xs[begin:end-1], xs[begin+1:end]) # (xβ,xβ), β¦, (xβββ,xβ)
sum(meth(f), pairs)
elseif method == "simpsons"
meth = f -> (lr -> begin l,r = lr; (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l) end)
xs = range(a, b, n+1)
pairs = zip(xs[begin:end-1], xs[begin+1:end]) # (xβ,xβ), β¦, (xβββ,xβ)
sum(meth(f), pairs)
elseif method == "midpoint"
xs = a:(b-a)/n:b
deltas = diff(xs)      # forms x2-x1, x3-x2, ..., xn-xn-1
cs = [(xs[i] + xs[i+1])/2 for i in 1:length(deltas)]
sum(f(cs[i]) * deltas[i] for i in 1:length(deltas))
end
end


but it works!!!

Maybe this is a bit more clear with how you can illustrate differences with various methods accompanying Riemann sums if performance isnβt the point:

function riemann(M, f, a, b, n)
xs = range(a,b,n)
riemann(M, f, xs[1:end-1], xs[2:end])
end
riemann(M, f, as, bs) = sum(M(f,a,b) * (b-a) for (a,b) in zip(as, bs))

Ms = (left = (f,a,b) -> f(a),
right= (f,a,b) -> f(b),
mid  = (f,a,b) -> f(a/2 + b/2),
mΜ    = (f,a,b) -> first(findmin(f, range(a,b,10))),
MΜ    = (f,a,b) -> first(findmax(f, range(a,b,10))),
trapezoid = (f,a,b) -> (f(a) + f(b))/2,
simpsons  = (f,a,b) -> (c = a/2 + b/2;(1/6) * (f(a) + 4*f(c) + f(b)))
)

ns = (5, 50, 500)

using DataFrames
d = DataFrame(n = Int[], M = Symbol[], value = Float64[])
for (m,M) in pairs(Ms)
for n in ns
v = riemann(M, sin, 0, pi, n)
push!(d, (n=n, M=m, value=v))
end
end

unstack(d, :n, :value)


Yielding with this f, a,b:

7Γ4 DataFrame
Row β M          5         50        500
β Symbol     Float64?  Float64?  Float64?
ββββββΌβββββββββββββββββββββββββββββββββββββββββ
1 β left        1.89612   1.99931   1.99999
2 β right       1.89612   1.99931   1.99999
3 β mid         2.05234   2.00034   2.0
4 β mΜ           1.11072   1.93523   1.9937
5 β MΜ           2.68152   2.06343   2.00629
6 β trapezoid   1.89612   1.99931   1.99999
7 β simpsons    2.00027   2.0       2.0

1 Like

I try to change ns = (5, 9, 17)

the ns actually represents the number of partitions + 1

when I see the solution manual for function 1/x it matches with the solution:

I have tried to add a new method in the Ms, basically I want to follow this formula:

T - \frac{(f'(b) - f'(a))h^{2}}{12}

with T is the result from the Trapezoidal, and h = \frac{b-a}{n}, n is the number of partition, and the interval is [a,b]

it seems the code still not produces the correct calculation, I shall try again:

...
modifiedtrapezoid = (f,a,b) -> (f(a) + f(b))/2 - (((f'(b)-f'(a))*h^2)/12),
...


Edited, this code works:

using CalculusWithJulia

# For partition of 8
ns = (9)

# Define the function, the partitions, and the interval
f(x) = x^4
a = 1
b = 3
n = 8
h = (b-a)/n
g = (((f'(b)-f'(a))*h^2)/12)
function riemann(M, f, a, b, n)
xs = range(a,b,n)
riemann(M, f, xs[1:end-1], xs[2:end])
end
riemann(M, f, as, bs) = sum(M(f,a,b) * (b-a) for (a,b) in zip(as, bs))

Ms = (left = (f,a,b) -> f(a),
right= (f,a,b) -> f(b),
midpoint  = (f,a,b) -> f(a/2 + b/2),
#mΜ    = (f,a,b) -> first(findmin(f, range(a,b,10))),
#MΜ    = (f,a,b) -> first(findmax(f, range(a,b,10))),
trapezoid = (f,a,b) -> (f(a) + f(b))/2,
modifiedtrapezoid = (f,a,b) -> (f(a) + f(b))/2 - g/2,
simpsons  = (f,a,b) -> (c = a/2 + b/2;(1/6) * (f(a) + 4*f(c) + f(b)))
)

using DataFrames
d = DataFrame(n = Int[], M = Symbol[], value = Float64[])
for (m,M) in pairs(Ms)
for n in ns
v = riemann(M, f, a, b, n)
push!(d, (n=n-1, M=m, value=v))
end
end

unstack(d, :n, :value)


The formula used to add the terms us: sum(M(f,a,b) * (b-a) for (a,b) in zip(as, bs)) has a factor h in the (b-a) term, so you only have one h in the formula, not h^2. Try

ct = (f,a,b)-> (c = a/2+b/2; (f(a)+f(b))/2 + 1/12 * (f'(b)-f'(a))*(b-a))

1 Like

Alright, thanks !