riemann-sum/riemann-sum.py

# Assignment Name: Integral Approximator
# Author: Rekai Nyangadzayi Musuka
# Description: Approximates Any Given Function
# Inputs: None
# Outputs: A Graph in Turtle with the Approximated and Theorized Area, highlighted by drawing trapezoids or squares.
# Version: 2019.0107.0

#LINK TO GRAPH (DESMOS): https://www.desmos.com/calculator/igcdhrepy3

import math
import turtle

CURSOR_SIZE = 0.1 # Makes the "Turtle" Invisible
DRAW_SPEED  = 0

#Cosmetic
FONT_FAMILY =  "Consolas"
NUM_FONT_SIZE = "8"
LABEL_FONT_SIZE = "10"
TEXT_SPACING = 15
FIGURES = 5

# Changes these to mess with the scale of the graph
SCALE_X     = 80 # Scale of the X Axis OG: 80
SCALE_Y     = 80 # Scale of the Y Axis OG: 80
UNIT_LENGTH = 400 # Determines the maximum domain and Range of the cartesian plane (in this case (-400, 400) in both X and Y)


# Best not to mess with these :)
MIN_VALUE   = -2 # The Domain of X is -2 <= x <= 2
MAX_VALUE   = 2
TICK_STEP   = 1 # Draw Ticks according to an interval of this variable
PLANE_INCREMENT = 1

# Change DRAW_STEP for a more or less accurately drawn function
DRAW_STEP   = 0.01 # What we Increment by when drawing the function

# Change to shorten or lengthen the height of the tick
TICK_HEIGHT = 10 # The Height of the Tick line

# Changes How Detailed the Area Approximation (Trapezoid or Rectangle) is (higher is better)
ITERATIONS = 100

# Determines whether the program uses the rectangle rule 
# or trapezoid rule (trapezoid is more exact I think?)

# TODO? Implement Simpson's Rule?
APPROX_TYPE = "trapezoid" # or square

# Initializing Turtle Variables
cursor = turtle.Turtle()
win = turtle.Screen()
cursor.shape()
cursor.hideturtle()
cursor.speed(DRAW_SPEED)

# Draws a Cursor Tick
def draw_tick(num, cursor, height, x_axis = False):
  
  # Function draws the numerical value of the x or y value at any given tick
  def draw_num(angle, distance, num):
    cursor.penup()
    cursor.left(angle)
    cursor.forward(distance)
    cursor.write(num, align="center", font=(FONT_FAMILY, NUM_FONT_SIZE, "normal"))
    cursor.backward(distance)
    cursor.left(-angle)
    cursor.pendown()

  if x_axis:
    draw_num(-90, (TICK_HEIGHT * 2), num)
  else:
    draw_num(90, (TICK_HEIGHT * 2), num)

  for angle in [90, -180]: 
    cursor.left(angle)
    cursor.forward(height / 2)
    cursor.backward(height / 2)
  cursor.left(90)
  

# Draws the Cartesian Plane
def create_plane():
  cursor.penup()

  cursor.goto(-UNIT_LENGTH, 0) # Goes to Fartheset left coordinate
  cursor.pendown()

  # Function will draw text nearby a given x value (used to draw axis labels).
  def draw_axis_label(angle, distance, text):
    cursor.penup()
    cursor.left(angle)
    cursor.forward(distance)
    cursor.write(text, align="center", font=(FONT_FAMILY, LABEL_FONT_SIZE, "normal"))
    cursor.backward(distance)
    cursor.left(-angle)
    cursor.pendown()

  # This draws the x line of the cartesian plane
  while (cursor.pos()[0]  <= UNIT_LENGTH): 
    if (cursor.pos()[0] == -UNIT_LENGTH): draw_axis_label(-90, TICK_HEIGHT * 4, "X Axis") # If x is the furthermost left value, draw the x axis label below
    if (cursor.pos()[0] % (TICK_STEP * SCALE_X) == 0): draw_tick(cursor.pos()[0], cursor, TICK_HEIGHT, True) # Draw ticks whenever the condition is met
    cursor.forward(PLANE_INCREMENT)
  
  cursor.penup()
  cursor.goto(0, -UNIT_LENGTH) # Goes to the farthest bottom coordinate
  cursor.pendown()
  cursor.left(90)
  
  # This draws the y line of the cartesian plane
  while (cursor.pos()[1]  <= UNIT_LENGTH):
    if (cursor.pos()[1] == -UNIT_LENGTH): draw_axis_label(90, TICK_HEIGHT * 4, "Y Axis") # If y is the furthermost bottom value, draw the y axis label there
    if (cursor.pos()[1] % (TICK_STEP * SCALE_Y) == 0): draw_tick(cursor.pos()[1], cursor, TICK_HEIGHT, True) # Draw ticks whenever the if condition is met
    cursor.forward(PLANE_INCREMENT)
  
  cursor.penup()

# Used for Calculus Portion of Program
def draw_rectangle(length, width):
  arr = [width, length, width, length]
  i = 0

  while (i < len(arr)):
    cursor.forward(arr[i])
    cursor.left(90)
    i += 1
  
  cursor.penup()
  cursor.forward(width)
  cursor.pendown()


# Used for Calculus Portion of Program
# Should Rewrite this sometime
def draw_trapezoid(a, b, h):
  # Drawing Trapezoid given Area
  width_of_triangle = b - a 

  c = math.sqrt(math.pow(h, 2) + math.pow(width_of_triangle, 2))

  angle = math.atan2(h, width_of_triangle) * (180 / math.pi)

  cursor.forward(h)
  cursor.left(90)
  cursor.forward(b)

  cursor.left(180 - angle)
  cursor.forward(c)
  cursor.left(angle)
  cursor.forward(a)
  cursor.left(90)

  cursor.penup()
  cursor.forward(h)
  cursor.pendown()


# The function that draws the upper half of the Heart
def f(x):
  return math.sqrt(1 - ((math.fabs(x) - 1) ** 2))

# The function that draws the Lower Half the the Heart
def g(x):
  return -3 * math.sqrt(1 - math.sqrt(math.fabs(x) / 2 ))

# Calls create_plane() which draws and labels the cartesian plane.
create_plane()


# This will set x to -2 and this while loop will then draw f(x) at -2 to 2
# Iterates between MIN_VALUE and MAX_VALUE and draws the calculated y value given the x value
x = MIN_VALUE
while (x <= MAX_VALUE):
  cursor.goto(SCALE_X * x, SCALE_Y * f(x)) # f(x) is the top half of the heart
  cursor.pendown()
  x += DRAW_STEP

cursor.penup()

# # This will set x to -2 and loop through g(x) until it reaches 2
# # Iterates between MIN_VALUE and MAX_VALUE and draws the calculated y value given the x value
x = MIN_VALUE
while (x <= MAX_VALUE):
  cursor.goto(SCALE_X * x, SCALE_Y * g(x)) # g(x) is the bottom half of the heart
  cursor.pendown()
  x += DRAW_STEP


# This code draws the functions used onto the screen.
cursor.penup()
cursor.goto(UNIT_LENGTH / 2, UNIT_LENGTH / 2) # We go to the middle of quadrant I
cursor.right(180)
cursor.write("Functions:", align="center", font=(FONT_FAMILY, NUM_FONT_SIZE, "normal"))
cursor.forward(TEXT_SPACING)
cursor.write("Top: y = math.sqrt(1 - ((math.abs(x) - 1) ** 2))", align="center", font=(FONT_FAMILY, NUM_FONT_SIZE, "normal"))
cursor.forward(TEXT_SPACING)
cursor.write("Bottom: y= math.sqrt(1 - math.sqrt(math.abs(x) / 2 ))", align="center", font=(FONT_FAMILY, NUM_FONT_SIZE, "normal"))

# Intergral Stuff Goes Here 

cursor.left(90) # Settings Cursor Back to it's default position (heading right)


# Allows us to import a reasonably acccurate estimate as to what the intergral of a given function will be
# By accurate I mean more accurate than our estimate.
import scipy.integrate as integrate

# Calculates Percent Error (the one we learned in chemistry)
def percent_error(exper, accept):
  return (math.fabs(exper - accept) / accept) * 100

# Approximates the intergral of a function using the square rule.
def approx_integral_square(func, a, b, n):
  # n is how many rectangles you want.
  # area is length * width
  area = 0.0
  width = (b - a) / n
  i = 1

  while (i <= n):
    height = func(a + (i - 1) * width)

    draw_rectangle(height * SCALE_Y, width * SCALE_Y)

    area += width * height
    i += 1

  return area

# Approximates the integral of a function using the trapezoid rule.
def approx_integral_trapezoid(func, a, b, n): 
  # n is how many trapezoids you want
  width = (b - a) / n 
  res = func(a) + func(b)
  i = 1
  
  prev_height = func(a)

  while (i < n):
    height = func(a + (i * width))
    draw_trapezoid(prev_height * SCALE_Y, height * SCALE_Y, width * SCALE_Y)
    res += 2 * height 

    prev_height = height
    i += 1

  return (width / 2) * res


# Exact Results
upper_area = integrate.quad(lambda x: f(x), MIN_VALUE, MAX_VALUE)
lower_area = integrate.quad(lambda x: g(x), MIN_VALUE, MAX_VALUE)

# Don't really care for the margin of error at the moment
# Absolute Value of the area because we don't care about signed area at the moment.math.
upper_area = math.fabs(upper_area[0])
lower_area = math.fabs(lower_area[0])

upper_area_approx = 0
lower_area_approx = 0

if (APPROX_TYPE == "square"):
  cursor.goto(MIN_VALUE * SCALE_X, 0)
  upper_area_approx = math.fabs(approx_integral_square(lambda x: f(x), MIN_VALUE, MAX_VALUE, ITERATIONS))
  
  cursor.goto(MIN_VALUE * SCALE_X, 0)
  lower_area_approx = math.fabs(approx_integral_square(lambda x: g(x), MIN_VALUE, MAX_VALUE, ITERATIONS))
elif (APPROX_TYPE == "trapezoid"):
  cursor.goto(MIN_VALUE * SCALE_X, 0)
  upper_area_approx = math.fabs(approx_integral_trapezoid(lambda x: f(x), MIN_VALUE, MAX_VALUE, ITERATIONS))
  
  cursor.goto(MIN_VALUE * SCALE_X, 0)
  lower_area_approx = math.fabs(approx_integral_trapezoid(lambda x: g(x), MIN_VALUE, MAX_VALUE, ITERATIONS))


# Probably Want to Display the Results or smth here.

cursor.penup()
cursor.goto(-(UNIT_LENGTH / 2), UNIT_LENGTH / 2) # We go to the middle of quadrant II
cursor.right(90)
if (APPROX_TYPE == "square"):
  cursor.write("Using the Square Rule:", align="center", font=(FONT_FAMILY, NUM_FONT_SIZE, "normal"))
else:
  cursor.write("Using the Trapezoid Rule:", align="center", font=(FONT_FAMILY, NUM_FONT_SIZE, "normal"))

cursor.forward(TEXT_SPACING)
cursor.forward(TEXT_SPACING)

cursor.write(("Approximated Top Half Area: " + str(round(upper_area_approx, FIGURES)) + " | Theorized: " + str(round(upper_area, FIGURES))), align="center", font=(FONT_FAMILY, NUM_FONT_SIZE, "normal"))
cursor.forward(TEXT_SPACING)
cursor.write(("Percent Error: " + str(round(percent_error(upper_area_approx, upper_area), FIGURES)) + "%"), align="center", font=(FONT_FAMILY, NUM_FONT_SIZE, "normal"))

cursor.forward(TEXT_SPACING)
cursor.forward(TEXT_SPACING)
cursor.write(("Approximated Bottom Half Area: " + str(round(lower_area_approx, FIGURES)) + " | Theorized: " + str(round(lower_area, FIGURES))), align="center", font=(FONT_FAMILY, NUM_FONT_SIZE, "normal"))
cursor.forward(TEXT_SPACING)
cursor.write(("Percent Error: " + str(round(percent_error(lower_area_approx, lower_area), FIGURES)) + "%"), align="center", font=(FONT_FAMILY, NUM_FONT_SIZE, "normal"))

win.exitonclick()