Problem

Integration by midpoint rule: The idea of the Midpoint rule for integration is to divide the area under a curve $f(x)$ into $n$ equal-sized rectangles. The height of the rectangle is determined by the value of $f$ at the midpoint of the rectangle. The figure below illustrates the idea,

To implement the midpoint rule, one has to compute the area of each rectangle, sum them up, just as in the formula for the Midpoint rule,

\[\int^b_a f(x) dx \approx h\sum^{n-1}_{i=0} f(a + i\times h + 0.5h) ~,\]

where $h=(b-a)/n$ is the width of each rectangle. Implement this formula as a Python function with this interface integrate(evalIntegrand, lowerLimit, upperLimit, numBin = 1000) where numBin represents $n$ in the above equation. Test the integrator with the following example input mathematical function,

import numpy as np
integrate(np.sin,0,2*np.pi,100)
-1.743934249004316e-17  

Note that this number could be slightly different on your system depending on your Python version and computer architecture, however, the number should be of the same order as the number in the above. Now, suppose you want to implement a projectile motion described via the physical equation $v = v_0 -g\times t$ where $v_0$ represents the initial velocity of the projectile, $g$ is the gravity constant, and $t$ is time starting from zero. Integrating this function from $t=0$ to $t=1$ would give us the total distance traveled by the projectile during this time period. Now, create a class implementation of this projectile motion such that the following code works and compute the total distance traveled by this projectile in the first second of the projectile motion.

projectile = Projectile(initVelocity=10)
print( "v(time=1) =", projectile.getVelocity(time=1) )
print( "Total distance traveled after 1 second: {} meters".format(integrate(projectile,0,1)))
v(time=1) = 0.1899999999999995  
Total distance traveled after 1 second: 5.095000000000001 meters  

Solution

Python

Here is an implementation of the integrate() function,

def integrate(evalIntegrand, lowerLimit, upperLimit, numBin = 1000):
    integral = 0
    binWidth = (upperLimit - lowerLimit) / float(numBin)
    for i in range(1, numBin + 1):
        integral += evalIntegrand(lowerLimit - 0.5 * binWidth + i * binWidth)
    integral *= binWidth
    return integral

Here is an implementation of the Projectile() function,

class Projectile(): 
    gravityConstant = 9.81
    def __init__(self, initVelocity):
        self.initVelocity = initVelocity
    def __call__(self, time):
        return self.getVelocity(time)
    def getVelocity(self,time):
        return self.initVelocity - self.gravityConstant * time

Comments