# Area che sta tra il grafico di una funzione limitata e continua a tratti e l'asse delle
# x, per x che varia tra a e b, approssimando il grafico con una poligonale ad N lati.
# Qui il grafico è fatto per punti, così da rappresentare correttamente anche le funz. discontinue
# Se aggiungo il comando  grafico <- 0  viene svolto solo il calcolo (per N grande conviene).
grafico <- 1
area <- function(F,a,b,N) {
   n <- N+1; x <- seq(b,a,len=n); y <- F(x)
   # suddivido [b,a] in N parti di estremi x[1]..x[n], prendo gli output y[i]
   # aggiungo x[n+1]=a e x[n+2]=b con ordinata sull'asse x: y[n+1]=y[n+2]=0
   x <- c(x,a,b); y <- c(y,0,0); n <- n+2
   if(grafico == 1) {
   plot(F,a,b, n=5000, type="p", pch=".", asp=1)
   abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3)
   for(i in 1:n) { points(x[i],y[i],pch=20,col="red")
     lines(c(x[i],x[i+1]),c(y[i],y[i+1]),col="red",lty=2)}
   }
   A <- (y[n]+y[1])*(x[n]-x[1])
   for(i in 1:(n-1)) A <- A + (y[i]+y[i+1])*(x[i]-x[i+1])
   A/2  # ho calcolato la somma delle aree degli N trapezi}
}
#
  
# esempio:
g <- function(x) ifelse(x<=2, x^2, x)
area(g,0,4, 1)
# Ho provato con la poligonale ad 1 lato. Ora aumento il n. dei lati:
area(g,0,4, 10)
area(g,0,4, 1e2)
area(g,0,4, 1e3)
# intuisco che è 8.6 e rotti; evito di rifare il grafico
grafico <- 0
area(g,0,4, 1e4)
area(g,0,4, 1e5)
area(g,0,4, 1e6)
# intuisco che potrebbe essere 8.666…