An archived instance of discourse for discussion in undergraduate Real and Numerical Analysis.

A Taylor approximation

mark

Estimate $\sin(1)$ to within $10^{-8}$ of the correct value using Taylor's theorem. We are interested in

  • An expansion of the sine around zero
  • An expansion of the sine around $\pi/3$ for comparison
  • Justification that our estimate is close enough
  • Python code that performs the computaiton
hjoseph

Last night I wrote 2 methods that approximate $sin(1)$ to within $10^{-8}$. I originally had a variable $fnxe$ in both loops that was based on the function, but because the series was centered at 0 I had to use an index of 2 which caused it to make an error in... the error measurement... When I decided to increment the count by 2 instead of 1, I skipped over all of the terms that have $sin(x)$ in them to avoid the $sin(0)$ terms because they caused the loop to breakout. This meant when I was on the 12th term:
$$\frac{-cos(x)x^{12}}{12!}$$
I skipped over the 13th term because n is incremented by 2, so instead of checking the error on:
$$\frac{sin(\xi)x^{13}}{13!}$$
I was checking it on:
$$\frac{cos(\xi)x^{14}}{14!}$$
I removed the $fnxe$ from the error term and replaced it with a 1, because $sin(x) \leq 1$ for all $x$. Since my error function is no longer dependent on derivatives of $sin(x)$, it is always nonzero and I can use an increment of $n= n+1$. It now uses only 12 steps as it should. The second loop is practically identical to the first mechanically, but is centered at $\frac{\pi}{3}$.





 """     
 The First method approximates using sin(0) in 12 steps.

 """

from numpy import pi
from sympy.mpmath import *
import math as math

y = 0    
n = 1
error = 1
while error > 1e-8:
    
   fnx = diff(sin,0,n)
   y = y + (fnx*(1)**n)/(math.factorial(n))
   error = abs((1)**(n+1))/(math.factorial(n+1))
   n = n+1
    
y,n

"""
The Second method approximates using sin(pi/3) in 5 steps. This method uses the produces the 
smallest taylor polynomial possible to be consistent to the 10^-8 decimal.
"""

from numpy import pi
from sympy.mpmath import *
import math as math

y = 0    
n = 0
error = 1
while error > 1e-8:

    fnx = diff(sin,pi/3,n)
    y = y + (fnx*(1-(pi/3))**n)/(math.factorial(n))
    error = abs((1-(pi/3))**(n+1)/(math.factorial(n+1)))
    n = n+1

y,n

Though you will never have too much error using just a 1 in place of $f^n(\xi)$, it seems like you may end up being too precise in some cases and generating a larger polynomial than you need. Is this possible? I know if you choose $f^n(\xi)>1$ it can certainly happen for $sin$ and $cos$ functions. I could see using a function that can decide whether or not to turn $sin$ into $csc$ or $cos$ into $sec$ or using an offset of $\frac{\pi}{4}$ so you could choose practically any $\xi$ based on proximity to $x$.

mark

@hjoseph I think this looks great! Could you (or maybe someone else) type out the mathematics involved using MathJax?

djett

Can someone explain the numpy.diff function? One source says it's function is to "Calculate the n-th order discrete difference along given axis". Any other way to describe that?

mark

@djett numpy.diff does not compute a derivative, though it is related. It computes differences between consecutive terms in an array. Thus,

from numpy import diff
diff([1,3,8,10])
#Out: array([2, 5, 2])

Note that @hjoseph used sympy.mpmath.diff, which is a different thing. Later this semester, we'll discuss scipy.misc.derivative, which is a numerical differentiator.

XCNate

So we know that $$f(x) = \sum_{k = 0}^n \frac{f^{(k)}(c)}{k!}(x - c)^k + \frac{f^{(n+1)}(\xi)}{(n+1)!}(x-c)^{n+1}.$$ In the case of the expansion $\sin(1)$ centered around $0$ our remainder will be $$ \frac{\sin^{(n+1)}(\xi)}{(n + 1)!}(1)^{n + 1} = \frac{\sin^{(n+1)}(\xi)}{(n + 1)!}.$$ Note that $\sin^{(n + 1)}(\xi) \leq 1$. Hence in order to find the correct $n$ to stay below our error bound we solve the inequality $$\left|\frac{\sin^{(n+1)}(\xi)}{(n + 1)!}\right| \leq \left|\frac{1}{(n + 1)!}\right| \leq 10^{-8}$$

y = 1
n = 0
while abs(y) > 1e-8:
    y = 1/math.factorial(n + 1)
    n = n + 1

n
#out: 12

In this case $n = 12$. So to approximate $\sin(1)$ we have

while n <= 12:
    y = y + diff(sin, 0, n)/math.factorial(n)
    n = n + 1

y
#out: mpf('0.84147098464806802')

Similarly, in the case of the expansion $\sin(1)$ centered around $\frac{\pi}{3}$ our inequality will be $$\left|\frac{\sin^{(n+1)}(\xi)}{n + 1}(1-\frac{\pi}{3})^{n + 1}\right| \leq \left|\frac{(1-\frac{\pi}{3})^{n+1}}{(n + 1)!}\right| \leq 10^{-8}$$

y = 1
n = 0
while abs(y) > 1e-8:
    y = ((1-(pi/3))**(n + 1))/math.factorial(n + 1)
    n = n + 1

n
#out: 5

In this case $n = 5$ so we have

y = 0
n = 0
while n <= 5:
   y = y + diff(sin, pi/3, n)/math.factorial(n)*(1-(pi/3))**n
   n = n + 1

y
#out: mpf('0.84147098482114002')
mark

@XCNate Looks cool - but I'm not seeing the output. Makes it hard to decide. I often include textual output as a comment, like:

2+2
#Out: 4
XCNate

Thanks for the suggestion! I was wondering what the best way to communicate the output was. I went ahead and edited the post to show the results.