This repository has been archived on 2021-06-30. You can view files and clone it, but cannot push or open issues or pull requests.
riemann-sum/riemann-sum.py

299 lines
10 KiB
Python

# 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()