-- integration.lua: vários métodos para integração numérica,
-- e alguns testes (para serem rodados em modo interativo).
-- Eduardo Ochs, 2010oct18
--
-- http://angg.twu.net/LUA/integration.lua.html
-- http://angg.twu.net/LUA/integration.lua
--         (find-angg "LUA/integration.lua")
-- http://angg.twu.net/2010.2-C2.html
--           (find-TH "2010.2-C2")
-- (find-eea2ps-links "")


pi = math.pi

particao_obvia = function (n_, x_0_, x_n_, crp)
    n, x_0, x_n = n_, x_0_, x_n_
    x_i = function (i) return x_0 + (x_n - x_0) * (i or i)/n end
    Critical_points = crp
  end
particao_explicita = function (X_, crp)
    n, X = #X_, X_
    x_i = function (i) return X[i] end
    Critical_points = crp
  end

-- x_i = function (i) return x_0 + (x_n - x_0) * (i or i)/n end
-- x_i = function (i) return X[i] end
x_i_1  = function (i) return x_i(i - 1) end
x_i_12 = function (i) return (x_i(i - 1) + x_i(i))/2 end
Dx     = function (i) return x_i(i) - x_i_1(i) end

y_i    = function (i) return f(x_i(i)) end
y_i_1  = function (i) return f(x_i_1(i)) end
y_i_12 = function (i) return f(x_i_12(i)) end
Dy     = function (i) return y_i(i) - y_i_1(i) end

fz_left    = function (i) return y_i_1(i) end
fz_middle  = function (i) return y_i_12(i) end
fz_right   = function (i) return y_i(i) end
fz_trapeze = function (i) return (y_i(i - 1) + y_i(i))/2 end
fz_simpson = function (i) return 1/3 * fz_trapeze(i) + 2/3 * y_i_12(i) end

riemann_sum = function (fz)
    local total = 0
    for i=1,n do total = total + fz(i) * Dx(i) end
    return total
  end

critical_points = function (i)
    C = {x_i_1(i), x_i(i)}
    for _,x in ipairs(Critical_points) do
      if C[1] < x and x < C[2] then table.insert(C, x) end
    end
    return C
  end
fz_inf = function (i)
    local C = critical_points(i)
    local y = f(C[1])
    for _,x in ipairs(C) do y = math.min(y, f(x)) end
    return y
  end
fz_sup = function (i)
    local C = critical_points(i)
    local y = f(C[1])
    for _,x in ipairs(C) do y = math.max(y, f(x)) end
    return y
  end

print_riemann_sums = function ()
    printf("inf:     %16.12f\n", riemann_sum(fz_inf))
    printf("sup:     %16.12f\n", riemann_sum(fz_sup))
    printf("left:    %16.12f\n", riemann_sum(fz_left))
    printf("right:   %16.12f\n", riemann_sum(fz_right))
    printf("middle:  %16.12f\n", riemann_sum(fz_middle))
    printf("trapeze: %16.12f\n", riemann_sum(fz_trapeze))
    printf("simpson: %16.12f\n", riemann_sum(fz_simpson))
  end


--[[

* (eepitch-lua51)
* (eepitch-kill)
* (eepitch-lua51)
ee_dofile "~/LUA/integration.lua"

f = function (x) return 4 - x^2 end
-- / x=pi
-- |      f(x) dx = 4 * 4 * 2/3 = 10.66666...
-- / x=0
particao_explicita({[0]=-2, -1, 1, 2}, {0})
print_riemann_sums()

f = function (x) return math.sin(x) end
-- / x=pi
-- |      f(x) dx = 2
-- / x=0
particao_obvia( 5, 0, pi, {pi/2}); print_riemann_sums()
particao_obvia(10, 0, pi, {pi/2}); print_riemann_sums()
particao_obvia(20, 0, pi, {pi/2}); print_riemann_sums()

f = function (x)
    if x <= 1 then return -1 end
    if x >= 3 then return 0 end
    return 4 - 2 * x
  end
eps = 1/1000000
crps = {1-eps, 1+eps, 3-eps, 3+eps}
particao_obvia( 20, 0, 4, crps); print_riemann_sums()
particao_obvia(100, 0, 4, crps); print_riemann_sums()
particao_obvia(500, 0, 4, crps); print_riemann_sums()
particao_obvia(501, 0, 4, crps); print_riemann_sums()

--]]

-- Local Variables:
-- coding:  raw-text-unix
-- modes:   (fundamental-mode lua-mode)
-- End: