Calculus 2, Lab 2: Arc length and numerical integration#
We have two goals for this tutorial assignment:
Understand how to use the arc length formula
Learn how to manipulate arrays and set up numerical integration
The second goal is closely tied to the first: most of the integrals that result when the arc length formula is applied cannot be evaluated by hand. We will begin by introducing arc length.
Let \(f(x)\) be a function that is differentiable on a closed interval \([a, b]\). (For the meaning of differentiability on a closed interval, see Section 2.1 in the textbook.)
The arc length of the curve \(y=f(x)\), for \(a\leq x\leq b\), is given by
The argument for this formula is quite clever: we partition the interval \([a, b]\) using points
and approximate the curve by the polygonal path passing through the points \((x_i, f(x_i))\), for \(i=0,\ldots, n\).
The length of each segment is given by the distance formula:
Now, by the Mean Value Theorem, on each interval \([x_i,x_{i+1}]\), there exists some \(c_i\) such that
giving us \(\Delta L_i = \Delta x_i \sqrt{1+(f'(c_i))^2}\). Putting this into a Riemann sum, we get our formula.
The trouble is that for all but the most carefully chosen curves, the arc length formula results in an integral that is impossible to evaluate symbolically. (Have a look at the functions that appear in the exercises for Section 7.4!)
To begin, make sure to run the cell below.
import sympy as sy
sy.init_printing()
x=sy.symbols('x')
Problem 1#
The curve \(y=x^2\) is usually considered to be one of the simplest examples in calculus.
(a) Use SymPy to define the function you need to integrate in order to compute the arc length. That is, first define the function \(f(x)=x^2\), and then define \(s(x) = \sqrt{1+f'(x)^2}\). Display the result to confirm it worked.
(b) What trig substitution could you use to evaluate the integral of this function? You don’t have to do the substitution; just indicate what you would use. Trig functions can be typed as \sin(x), \cos(x), etc. Run the cell when you’re done to format the text.
(c) Use SymPy to calculate the length of the curve \(y=x^2\) from \(x=0\) to \(x=2\).
# The command to integrate f from a to b is sy.integrate(f,(x,a,b))
# You're computing arc length. Be sure to use the correct function in your integral!
Problem 2#
Consider problem 7.4.19 in the textbook, which advises you to set up, but not evaluate, the integral to compute the length of the curve \(y=1/x\), from \(x=1\) to \(x=2\).
(a) Use SymPy to define the function you will need to integrate, according to the arc length formula.
(b) Try to evaluate the integral as you did in Problem 1. Attempt to compute both the definite and indefinite integrals.
Does your answer not look particularly familiar? You can get a decimal approximation by adding .evalf() to the end of your definite integral:
Problem 3#
For this problem we will use the implementation of Simpson’s Rule given by Patrick Walls of UBC. Set things up as indicated, and then (approximately) evaluate the integral from Problem 2 using Simpson’s Rule.
import numpy as np
Here is Patrick’s code for a Simpson’s Rule function:
def simps(f,a,b,N=50):
'''Approximate the integral of f(x) from a to b by Simpson's rule.
Simpson's rule approximates the integral \int_a^b f(x) dx by the sum:
(dx/3) \sum_{k=1}^{N/2} (f(x_{2i-2} + 4f(x_{2i-1}) + f(x_{2i}))
where x_i = a + i*dx and dx = (b - a)/N.
Parameters
----------
f : function
Vectorized function of a single variable
a , b : numbers
Interval of integration [a,b]
N : (even) integer
Number of subintervals of [a,b]
Returns
-------
float
Approximation of the integral of f(x) from a to b using
Simpson's rule with N subintervals of equal length.
Examples
--------
>>> simps(lambda x : 3*x**2,0,1,10)
1.0
'''
if N % 2 == 1:
raise ValueError("N must be an even integer.")
dx = (b-a)/N
x = np.linspace(a,b,N+1)
y = f(x)
S = dx/3 * np.sum(y[0:-1:2] + 4*y[1::2] + y[2::2])
return S
Note that most of the above consists of comments: the actual code is only a few lines long!
You can run this by simply entering simps(f,a,b,N), where a,b,N are values you choose, and f is the function you want to integrate. For example, you might run simps(f,0,2,50) to integrate a function f from 0 to 2 using 50 subintervals.
Caution: this will not work if f is a SymPy function! You can either redefine the function f you want to integrate (and this needs to be the arc length expression, not the function whose length you’re computing), or convert the function from Problem 2 into a NumPy function. If you want to know how to do this, see the numerical integration example I prepared.
Compute the integral for the arc length of \(y=1/x\) from 1 to 2 using Simpson’s rule, with 10, 50, and 100 subintervals.
# Use SymPy to calculate the arc length function s(x)
# Define s(x) as a NumPy function
# Simpson's with 10 subintervals
# Simpson's with 50 subintervals
# Simpson's with 100 subintervals
How well do your results compare to what you found in Problem 2?
# Compute the difference of each answer above with the evalf() value from Problem 2
Problem 4#
Finally, let’s compute the length of a half circle. Consider the function \(f(x) = \sqrt{4-x^2}\) on the interval \([-2,2]\). We know that the circumference of a semi-circle of radius \(R\) is \(\pi R\), so without any calculus, the answer to the problem should be \(2\pi\).
We could try to compute this using Simpson’s rule, but there’s a problem: \(f'(x)\) is undefined at the endpoints, and Simpson’s rule uses the endpoints!
(Feel free to try it.)
Instead, we could try a Riemann sum using midpoints, to stay away from the endpoints. How do we set up a Riemann sum? First, we need to understand how to work with arrays in Numpy. Let’s start by defining an array, and since we’ll be needing it anyway, let’s make it a uniform partition.
P = np.linspace(0,10,11)
P
You should see an array consisting of the numbers from 0 to 10. Python has a “colon syntax” for accessing parts of an array. To see the full array P, you can enter P[:]. Try it:
If you want all the entries after a certain point in the index, you can place a colon to the right of the index value. For example, P[3:] will give everything after (and including) the third index position. Try it:
P[:3]
Remember that a Python index starts at 0, so position 3 is actually the 4th entry!
To get everything before a certain position, put the colon on the left. See what happens if you enter P[:4]:
Finally, you can also start from the end instead of the beginning, using negative index values. See what happens if you input P[:-3] and P[-4:]:
Why do we need this? Remember that to do a Riemann sum with left endpoints, we need all points in our partition except the last one. The set of left endpoints is given by:
L=P[:-1]
L
R=P[1:]
R
To get midpoints, we need to average the left and right endpoints. See what happens when you run the next cell:
M = (L+R)/2
M
Now we can set up a midpoint Riemann sum. We will try to imitate Patrick’s code for Simpsons rule.
def midpoint(f,a,b,N=50):
'''Approximate the integral of f(x) from a to b using the midpoint rule.
Parameters
----------
f : function
Vectorized function of a single variable
a , b : numbers
Interval of integration [a,b]
N : (even) integer
Number of subintervals of [a,b]
Returns
-------
float
Approximation of the integral of f(x) from a to b using
midpoint rule with N subintervals of equal length.
Examples
--------
>>> midpoint(lambda x : 3*x**2,0,1,10)
1.0
'''
dx = (b-a)/N
x = np.linspace(a,b,N+1)
mid = (x[:-1]+x[1:])/2
y = f(mid)
S = np.sum(y) * dx
return S
Now we can finally return to the problem we needed to solve. First, use SymPy to calculate \(\sqrt{1+f'(x)^2}\) for \(f(x)=\sqrt{4-x^2}\).
Next, either use the “lambdify” process to turn this into a NumPy function, or simply define the result as a NumPy function by hand.
Finally, use the midpoint sum function above to approximate the integral from \(-2\) to \(2\) using 10, 100, and 1000 rectangles. How well do your answers compare to the expected result of \(2\pi\)?