Week 11: Analyzing the Error Using Symbolic Method (SymPy)#
Laboratory 4e
Last updated July 26, 2022
00. Content #
Mathematics
Error analysis
Limits
Integrals
Derivatives
Programming Skills
Number type: float
Compute limits, integrals, and derivatives using SymPy
Embedded Systems
N/A
0. Required Hardware #
N/A
Copy the file `student_info.py` into the folder you are working in, then run this cell to put your name and email into this document.
from student_info import show_info
show_info()
1. Floating Point Numbers #
Unfortunately, it’s not possible to represent every number exactly in a computer. For numbers which are not whole numbers, Python only keeps a certain number of digits. This way of storing numbers is called floating point, and numbers stored this way are called floating point numbers. In Python, the type for floating point numbers is float.
Floating point numbers have some quirks you need to be aware of. For example, we know that 0.1+0.2=0.3, but look what Python tells us:
0.1+0.2
It is slightly wrong. What is going on here? Floating point numbers are stored in binary, and in binary, 0.1 and 0.2 are infinite repeating decimals. In binary, we would write the number one tenth as 0.0001100110011… and so on forever. We don’t have an infinite amount of memory available on our computers, so Python only stores 53 binary digits, which introduces roundoff error. The relative error is at most \(2^{-53} = 1.11\times 10^{-16}\), so numbers are not reliable beyond sixteen places. Sure enough, in the example above, we got an error sixteen places after the 3.
2**-53
Exercise 1 (5 pts)#
What happens if you multiply 0.1 and 0.2 by 10, add them, and then divide by 10? Mathematically, \(0.1+0.2 = \frac{0.1(10)+0.2(10)}{10}\), so do you still get a slightly wrong answer with floating point numbers?
Write Answers for Exercise 1 Below
# Type out the equation to see its output
2. Symbolic Computation #
Most of the time, working with numbers in a computer, close enough is good enough. We use floating point numbers, and they almost never give us any trouble. There are times, though, that you want the exact answer. That’s when you need symbolic computation. For this lab, instead of Numpy (short for “numerical Python”) we will use Sympy (“symbolic Python”). Let’s compare what happens when we ask each of them for a square root:
import numpy as np
import sympy
display(np.sqrt(12))
display(sympy.sqrt(12))
Numpy gets us the floating point representation of \(\sqrt{12}\) in no time flat. Sympy instead keeps an exact answer, though it kindly simplifies it to \(2\sqrt{3}\). Sympy is more than capable of providing a decimal approximation, to whatever level of precision you like:
sympy.sqrt(12).evalf(100)
Of course, you and I will probably never need 100 digits of precision for anything, but if we ever did need it, it’s there. The real power of Sympy is that you don’t need to ever work with a decimal representation. For example, let’s compare how Numpy and Sympy handle this procedure:
for x in [np.sqrt(12), sympy.sqrt(12)]:
display((x+1)*(x-1))
We are already seeing a little bit of floating point error creep in with Numpy due to the cancellation in the x-1 term. Sympy gets even more attractive when we use its ability to simplify:
((x+1)*(x-1)).simplify()
It got us the exact answer, no floating point error at all.
Formulas#
As you saw in lab 4d, we can also use Sympy to manipulate formulas. This can make it way easier to work with big expressions and prevent errors. It’s usually pretty simple to do, too. All we need is to create a variable using Sympy. There are two ways to do this. For common variables, you can use the sympy.abc module:
from sympy.abc import x, y, r, theta, phi
display(x, y, r, theta, phi)
For more complicated or uncommon symbols, you can use the function sympy.symbols to specify them using a bit of Latex.
e_r = sympy.symbols(r'\varepsilon_r')
e_theta = sympy.symbols(r'\varepsilon_\theta')
display(e_r, e_theta)
Notice we used a “raw string” with the r'...' syntax, so that we don’t have to worry about Python capturing any backslashes (backslash is also an escape character in Python, and it does different things from Latex).
Now that we have the variables, we can just start doing stuff with them.
fraction = (x**2 - y**2) / (x-y)
display(fraction, fraction.simplify())
When you do any operation with a Sympy object, the result becomes another Sympy object. That lets us use normal Python syntax without anything extra.
# x is a Sympy object, so x*2 is as well
display(x*2)
# then when we divide that by 3, we get another Sympy object
display((x*2)/3)
# and because of the order of operations, that's the same as if we dropped the parentheses
display(x*2/3)
You have to watch out, though: if you combine things which are not Sympy objects, Python will do its normal thing with them. That’s why this leads to a floating point in your answer:
x*(2/3)
Usually you can get around this by changing some parentheses. If necessary, you can also use sympy.S to convert any number into a Sympy number.
sympy.S(2)/3
Exercise 2 (10 pts)#
The expression below (((22/7)**2 - sympy.pi**2)**(1/2)) has multiple floating point numbers which snuck into the computation. Rewrite it so that there are no floating points.
Write Answers for Exercise 2 Below
((22/7)**2 - sympy.pi**2)**(1/2)
3. A Practical Computation #
Now that you know the basics of Sympy, let’s work together to use it to answer a potentially-tricky question which came up in lab 9. We have an expression for the error in our measurement, and we want to know how we should do the experiment to minimize the error. We can start with the formula.
y = r/(sympy.cot(theta)-sympy.cot(phi))
y
We can ask it for derivatives using builtin functions (I made a few errors when doing this one by hand earlier. It definitely makes sense to use a computer to check the answers here).
y.diff(theta)
y.diff(phi)
Then make an expression for the combined error.
y_error = e_r * y.diff(r) + e_theta * (y.diff(theta) - y.diff(phi))
y_error
That looks a bit messy. Let’s have the computer simplify it before we go on.
y_error = y_error.simplify()
y_error
Here is something which stands out about this result: the formula for y appears in our answer for the error. Let’s divide it out to get the relative error.
relative_error = (y_error/y).simplify()
relative_error
A bit nicer. Now, here is where we can do something that would be a real pain by hand but is easy by computer. I don’t like that the answer depends on \(\theta\) and \(\phi\), since I can’t control those. I would rather have an answer which depends on \(y\) (which is fixed) and \(z\) and \(r\) which I can control. Let’s do it:
from sympy.abc import y, z
relative_error = relative_error.replace(sympy.cot(phi), z/y)
relative_error
Exercise 3 (10 pts)#
Use a similar strategy to replace \(\cot\theta\) with an expression involving only \(y\), \(z\), and \(r\). After simplifying, you should get
Write Answers for Exercise 3 Below
4. Analyzing the Result #
Now that we have a (remarkably tidy!) expression for the relative error, we should revisit the things we can control. We could reduce \(\varepsilon_\theta\) by aligning and steadying the inclinometer more carefully. We could reduce \(\varepsilon_r\) by measuring our rope in a more accurate way, or using a more precise tool for measuring horizontal distance. What about \(r\) and \(z\)? How does \(z\) affect the answer?
Exercise 4 (10 pts)#
In two sentences, describe qualitatively how the relative error in the height depends on the parameter \(z\). Are bigger or smaller values better, as far as reducing error goes?
Write Answers for Exercise 4 Below
\(r\) dependence#
The dependence on \(r\) is a bit trickier. Let’s ask Sympy to tell us what happens as \(r\) approaches infinity.
relative_error.limit(r, sympy.oo)
The answer depends on the sign of \(\varepsilon_\theta\) and \(y\), but we know both are positive so this indicates that the error would grow to infinity if we used an infinitely long rope. Not good!
What about using a short rope?
relative_error.limit(r, 0)
Same thing! if we use a very short rope, we get infinite error. By the open interval method, we know there has to be an ideal rope length in between. Let’s find it by setting the derivative to zero. We can take the derivative:
derivative = relative_error.diff(r)
derivative
and the function sympy.solve will solve the equation derivative == 0 for \(r\).
solutions = sympy.solve(derivative, r)
solutions
The first one is negative, and our rope has a positive length. Let’s take the other one:
ideal_length = solutions[1].simplify()
ideal_length
Exercise 5 (20 pts)#
Use the ideal_length.subs and one of your measurements from lab 8 to figure out the ideal length of rope. You may need to compute \(z\), but remember that it is just \(y\cot\phi\). How does it compare to the length of rope you actually used? How would you measure it differently, if you wanted to go back and get a more accurate measurement?
Write Answers for Exercise 5 Below
5. Identifying Limiting Factors #
By putting our ideal length back into the expression for the error, we can get see how good our error can get. This will tell us where we should really spend our efforts if we want to get a really accurate answer.
best_error = relative_error.subs(r, ideal_length).simplify()
best_error
Some times, Sympy does not find the best answer right away. We can give it a kick in the right direction, though. Since we have a square root inside that square, it would probably get simpler if we expanded the square, then combined like terms. Let’s ask for that explicitly:
best_error = best_error.expand().simplify()
best_error
Very nice! The best error we can hope for, picking an ideal length of rope (after a bit more tidying) is
This lets us definitively say: the error depends really strongly on \(\varepsilon_\theta\). If we want a low relative error, we have to work to get that number down, since it multiplies into every other term. It only makes sense to reduce \(\varepsilon_r\) if the term \(\frac{\varepsilon_{r} }{\varepsilon_{\theta} y}\) is close to 2, otherwise it will be completely dwarfed.
Exercise 6 (10 pts)#
\(\varepsilon_\theta\) also appears in the denominator. Does that mean we should try and keep it from getting too small? Use Sympy to find \(\displaystyle\lim_{\varepsilon_\theta \to 0}\)best_error and use the result of that computation to justify a decision.
Write Answers for Exercise 6 Below
6. Do Your Homework Automatically #
In the example above we used lots of useful Sympy tools – substitution, simplification, differentiation, and taking limits were all essential. There are way more things Sympy can do besides this, which you can read about in its documentation. What you know already is enough to do some things which could really help you out when doing you calculus homework. For example, look at this question which was on the third midterm for MA161 in 2021:
With Sympy we check our work effortlessly!
function = (sympy.exp(3*x**2) -1)/(sympy.cos(x)-1)
display(function)
function.limit(x, 0)
There are tools to compute Taylor series as well, which might be helpful to check your work studying for the final in MA162. Take this problem, which was on the MA162 final in the fall of 2021:
Find the Taylor Series for \(f (x) = \sin(πx)\) centered at \(a = 1\).
sympy.sin(sympy.pi*x).series(x, 1)
Or you could do this integral, from the same exam.
sympy.integrate(sympy.exp(3*x)*sympy.sin(x), x).simplify()
Hopefully this has convinced you that Sympy can be a handy calculator to have around.
Exercise 7 (10 pts)#
Pick a problem from a previous calculus exam. Copy it here using Latex, and solve it using Sympy.