Integration
Finding the functions from derivatives is called integration.
""" Integrals / Integration
By antidifferention or integration we can get the original function
from the derived function.
Integration can be used to find areas, volumes, central points
and many other useful things.
Function: y = 3x + 2
Derivative: d = 3
Integral: I = 3x
y' = ax then y = (ax^2)/2 + C
y = ax + b then y = (ax^2)/2 + bx + C
"""
import numpy as np
import matplotlib.pyplot as plt
from sympy import *
import scipy.integrate as spi
# -----------------------------------------------------
x = Symbol('x')
t = Symbol('t')
f = 2*x
d = f.diff(x)
d_integration = integrate(d, x)
assert d_integration == f
# -----------------------------------------------------
# Function to be integrated
def func(x):
return 2*x
# Integration with scipy quad() - area under the f(x)
xa, xb = 0, 2
func_integral, err = spi.quad(func, xa, xb) # lower & upper limits
A = func_integral
# -----------------------------------------------------
fig, ax = plt.subplots()
x = np.linspace(0, 3, 100)
y = func(x)
ax.plot(x, y, label="f(x) = 2x")
ax.fill_between(x, y, 0, where=(x >= xa) & (x <= xb),
color="gray", alpha=0.5, label=f'Integral (Area) = %s' % A)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.legend()
plt.show()
print("Function f(x) =", f)
print("Derivative f' =", d)
print("Integral I =", d_integration)
print('Area A =', func_integral)
"""
Function f = 2*x
Derivative f' = 2
Integral I = 2*x
Area A = 4.0
"""
Definite integral
A definite integral has start and end values.
""" Integrals / Definite integral
A definite integral has start and end values [a, b]
The symbol for integral is a stylish S from sum, or summing slices
"""
import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
# ---------------------------------------------------------------
# Function to be integrated
def f(x):
return x**2
def integral_approximation_A(Y, a, b): # numpy
approx = (b-a) * np.mean(Y)
return approx
def integral_approximation_B(f, a, b): #scipy
approx, error = spi.quad(f, a, b)
return approx
# Integral boudaries (lower and upper limits)
a = 0
b = 1
# Function values
range = np.arange(a, b + .0001, .0001)
Y = f(range)
# Integral approximation
approx_A = integral_approximation_A(Y, a, b)
approx_B = integral_approximation_B(f, a, b)
# ---------------------------------------------------------------
print(approx_A) # 0.33335
print(approx_B) # 0.33333333333333337
# Plotting the function
x = np.linspace(0, 1, 100)
y = f(x)
fig, ax = plt.subplots()
ax.plot(x, y, label="f(x) = x^2")
ax.fill_between(x, y, 0, where=(x >= 0) & (x <= 1), color="gray", alpha=0.5,
label=f'I (Area) = %s' %round(approx_B, 4))
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.legend()
plt.show()
Gaussian function
Integrating the standard Gaussian distribution.
""" Integrals / Gausssian distribution
Gaussian distribution is also called normal distribution.
We calculate integral for gaussian functions using scipy approximation.
The definitive integral of a function gives the total area
under the curve function over an interval (a, b).
PDF means probability distribution function.
"""
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import scipy.stats as stats
# ---------------------------------------------------------------------
# Gaussian function
def gaussian(x, mu, sigma):
return np.exp((-1*(x - mu)**2) / (2 * sigma**2))
# Define mu and sigma
mu, sigma = 0, 1
# Axis for plot
X = np.arange(-5, 5, 0.001)
Y = stats.norm.pdf(X, mu, sigma) # calculates PDF for normal distribution
# Define integral boundaries
a = -3
b = 3
# Integral approximation (area)
approx, error = integrate.quad(
lambda x: gaussian(x, mu, sigma), a, b
)
# ---------------------------------------------------------------------
print("Gaussian distribution Area =", approx)
plt.plot(X, Y, label="Gaussian distribution")
plt.fill_between(X, Y, 0, where=(X >= a) & (X <= b), color='gray', alpha=0.5,
label=f'I = %s' %round(approx, 4))
plt.legend()
plt.show()
"""
Gaussian distribution Area = 2.4998608894830947
"""
Falling object
Use antidifferentiation to get distance traveled from acceleration.
""" Object Free fall
What is the time in which an object dropped from a 400ft point
reaches earth (in the absence of air)
t = ?
Galileo: An object is falling to Earth with the same acceleration:
a = 32 ft/s^2
a = 9.8 m/s^2
Acceleration is the instantaneus rate of change of speed / time
v' = 32
v = 32t + C
v = 32t # object is dropped, zero speed at start
Instantaneus speed is the rate of change of distance / time
s' = 32t
s = 16t^2 + C
s = 16t^2
To answer the initial question, 400ft falling object
400 = 16t^2
t = [-5, 5] # both -5 and 5 are correct
t = 5 # we choose 5, because we have downward fall
By antidifferentiation we can proceed from acceleration to speed,
and from acceleration to the distance traveled.
a = 32
v = a'
s = a''
Distance traveled in 5 seconds is 400
t = 5
s = 400
"""
import numpy as np
from sympy import *
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# ------------------------------------------------------------
t = Symbol('t')
a = 32
v = integrate(a, t)
s = integrate(v, t)
# Distance traveled in 5 seconds
d = s.subs(t, 5)
# Time for 400ft fall
s = s - 400 # s = 16*t**2 - 400
time = solve(s, t) # resolve equation that makes s = 0
# ------------------------------------------------------------
# Print results
print('Acceleration =', a) # Acceleration = 32
print('Speed =', v) # Speed = 32*t
print('Distance =', s) # Distance = 16*t**2
print("Distance traveled in 5 seconds =", d) # 400
print("Time for 400ft fall =", time[1]) # 5
# Plotting
t = np.linspace(0, 6)
s = 16*t**2
fig, ax = plt.subplots()
ax.plot(t, s, label="s(t) = 16t^2")
ax.set_xlabel("t")
ax.set_ylabel("s(t)")
ax.legend()
plt.scatter(5, 400, label="s(5) = 400")
plt.plot((5, 5), (0, 400), linestyle='--')
plt.plot((0, 5), (400, 400), linestyle='--')
plt.show()
# Animation
def update(frame):
t = np.linspace(0, frame/10)
s = 16 * t**2
ax.clear()
ax.plot(t, s)
ax.set_xlabel("t")
ax.set_ylabel("s(t)")
ax.set_ylim(400, 0)
ax.set_xlim(0, 5.2)
ax.set_title("s(t) = 16t^2")
y = 16*(frame/10)**2
ax.scatter(5, y, label="s(t) = 400 - 16t^2")
ax.plot((frame/10, frame/10), (5, y), linestyle='--')
ax.plot((5, frame/10), (y, y), linestyle='--')
fig, ax = plt.subplots()
ani = FuncAnimation(fig, update, frames=np.arange(10, 51, 1), repeat=True)
plt.show()
if 1 == 0: # save image
ani.save('1427_falling_ball.gif', writer='imagemagick', fps=10)
"""
Acceleration = 32
Speed = 32*t
Distance = 16*t**2 - 400
Distance traveled in 5 seconds = 400
Time for 400ft fall = 5
"""
Throwing down
Use integration to get time for 400ft thrown fall. When the object is trown down it has an starting velocity.
""" Object thrown downward
When the object leaves the hand it already has some velocity
Suppose, for example, that the velocity is 100 ft/sec
a = 32
v = 32t + C; 100 = 32 * 0 + C; C = 100
v = 32t + 100
s' = 16t^2 + 100t + C
Let us agree that the distance is measured from the release point
t = 0; s = 0; C = 0
s = 16t^2 + 100t
Distance traveled in 5 seconds is 400
t = 5
s = 900
"""
import numpy as np
from sympy import *
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# ------------------------------------------------------------
t = Symbol('t')
a = 32
v = integrate(a, t) + 100 # Look Hee
s = integrate(v, t)
# Distance after 5 sec
distance1 = s.subs(t, 5)
# Distance after 2.77 sec
distance2 = s.subs(t, 2.77)
# Time for 400ft thrown fall
s = s - 400 # s = 16*t**2 + 100t
time = solve(s, t) # resolve equation that makes s = 0
# ------------------------------------------------------------
# Print results
print('Acceleration =', a) # 32
print('Speed =', v) # 32*t + 100
print('Distance =', s) # 16*t**2 + 100*t
print("Distance traveled after 5 sec =", distance1) # 900
print("Distance traveled after 2.77 sec =", distance2) # 400
print("Time for 400ft thrown fall =", round(time[0], 2)) # 3
# ------------------------------------------------------------
# Plotting
t = np.linspace(0, 5)
s = 16*t**2 + 100*t
fig, ax = plt.subplots()
ax.plot(t, s, label="s(t) = 16t^2 + 100t")
ax.set_xlabel("t")
ax.set_ylabel("s(t)")
ax.legend()
plt.scatter(2.77, 400, label="s(.277) = 400")
plt.plot((2.77, 2.77), (0, 400), linestyle='--')
plt.plot((0, 2.77), (400, 400), linestyle='--')
plt.show()
# Animation
def update(frame):
t = np.linspace(0, frame/10)
s = 16 * t**2 + 100*t
ax.clear()
ax.plot(t, s)
ax.set_xlabel("t")
ax.set_ylabel("s(t)")
ax.set_ylim(400, 0)
ax.set_xlim(0, 5.2)
ax.set_title("s(t) = 16t^2 + 100t")
y = 16*(frame/10)**2 + 100*(frame/10)
ax.scatter(2.77, y, label="s(t) = 400 - (16t^2 + 100t)")
ax.plot((frame/10, frame/10), (2.77, y), linestyle='--')
ax.plot((2.77, frame/10), (y, y), linestyle='--')
fig, ax = plt.subplots()
ani = FuncAnimation(fig, update, frames=np.arange(10, 29, 1), repeat=False)
plt.show()
if 1 == 0: # save image
ani.save('1427_throw_ball.gif', writer='imagemagick', fps=10)
"""
Acceleration = 32
Speed = 32*t + 100
Distance = 16*t**2 + 100*t - 400
Distance traveled after 5 sec = 900
Distance traveled after 2.77 sec = 399.766400000000
Time for 400ft thrown fall = 2.77
"""
Last update: 370 days ago
