For CalculusWithJulia, is there Riemann' Midpoint method?

General Usage
1

I read this:

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
2

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
3

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?

4

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.

5

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
6

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!!!

ch1-calculus-5thedefiniteintegral-59

7

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
8

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:

Capture d’écran_2022-12-19_11-16-06
Capture d’écran_2022-12-19_11-16-06748Γ—210 26.5 KB

9

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)
10

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
11

Alright, thanks !