# Estimating the area under the curve x^2 using random points

# Given function
parabola <- function(x)
{
    return(x^2)
}
plot(parabola)

# Definite integral over the interval [a,b]
areaintegrparab <- function(a,b)
{
    return(b^3/3-a^3/3)
}

# Simulation approach to the calculation of the area
# a is the lower bound of [a,b]
# b is the upper bound of [a,b]
# c is the number of random points to be used
# e is the maximum allowed error (in percentage points)
areamt <- function(a,b,c,e)
{
    randx<- runif(c,a,b)
    randy<- runif(c,min(parabola(a),parabola(b)),max(parabola(a),parabola(b)))
    c <- c(0)
    h = 1
    n = 0
    m = 1
    o = 1
    vecy = c(0)
    vecx = c(0)
    for(number in randy)
    {
        if(number< parabola(randx[h]))
        {
            vecy[m]<-number
            vecx[o]<-randx[h]
            n = n + 1
            h = h + 1
            m = m + 1
            o = o + 1
        }else
        {
            h = h + 1
        }
    }

    # Points ratio
    pr=n/length(randx)

    # Estimated area
    mca = n/length(randx)*abs(b-a)*abs(max(parabola(a),parabola(b))-min(parabola(a),parabola(b)))

    # Exact area
    aint = areaintegrparab(a,b)

    # Error
    change = (mca-aint)/aint*100

    # Vector containing the summary of the entire process
    data = c(pr,mca,aint,change)

    # Let's put some labels
    if(change>=0)
    {
        namechange = "% Overestimated"
    }else
    {
        namechange = "% Underestimated"
    }
    names(data)=c("Points ratio","Estimated area","Real area",namechange)

    # Check if error is compatible with the maximum error allowed
    # If TRUE, data is printed and plotted, else the process is reapeted with more points
    if(abs(change) < e)
    {
        print(data)#Estimating the area under the curve x^2 using random points

# Given function
parabola <- function(x)
{
    return(x^2)
}
plot(parabola)

# Definite integral over the interval [a,b]
areaintegrparab <- function(a,b)
{
    return(b^3/3-a^3/3)
}

# Simulation approach to the calculation of the area
# a is the lower bound of [a,b]
# b is the upper bound of [a,b]
# c is the number of random points to be used
# e is the maximum allowed error (in percentage points)
areamt <- function(a,b,c,e)
{
    randx<- runif(c,a,b)
    randy<- runif(c,min(parabola(a),parabola(b)),max(parabola(a),parabola(b)))
    c <- c(0)
    h = 1
    n = 0
    m = 1
    o = 1
    vecy = c(0)
    vecx = c(0)
    for(number in randy)
    {
        if(number< parabola(randx[h]))
        {
            vecy[m]<-number
            vecx[o]<-randx[h]
            n = n + 1
            h = h + 1
            m = m + 1
            o = o + 1
        }else
        {
            h = h + 1
        }
    }

    # Points ratio
    pr=n/length(randx)

    # Estimated area
    mca = n/length(randx)*abs(b-a)*abs(max(parabola(a),parabola(b))-min(parabola(a),parabola(b)))

    # Exact area
    aint = areaintegrparab(a,b)

    # Error
    change = (mca-aint)/aint*100

    # Vector containing the summary of the entire process
    data = c(pr,mca,aint,change)

    # Let's put some labels
    if(change>=0)
    {
        namechange = "% Overestimated"
    }else
    {
        namechange = "% Underestimated"
    }

    names(data)=c("Points ratio","Estimated area","Real area",namechange)

    # Check if error is compatible with the maximum error allowed
    # If TRUE, data is printed and plotted, else the process is reapeted with more points
    if(abs(change) < e)
    {
        print(data)
        plot(parabola,xlim=c(a,b),ylim=c(parabola(a),parabola(b)),col="red",main="Estimating the area under the curve")
        points(vecx,vecy,col="green")
        return(n/length(randx)*abs(b-a)*abs(max(parabola(a),parabola(b))-min(parabola(a),parabola(b))))
    }else
    {
        return(areamt(a,b,c+1000,e))
    }
}

# Calculate the area under the curve between 0 and 2 with 30000 points
# and a maximum allowed error of 2%
areamt(0,3,30000,2)