🌱 Taylor series

we're back with more bad sines: order=5
a=0.000
reshape x as triangle wave

this time we've got controools! you're looking at a "taylor series" approximating sin(x). loosely speaking, the order is how many steps to take in the approximation (or, the number of iterations in the for-loop, minus 1), and a is the "anchor point", around which the approximation is most accurate (notice how the error plot drops around it). the checkbox enables the triangle-wave symmetry-hack from last post.

before we continue you may have a complaint about the code:

isn't it weird that the approximation of sin(x) itself contains sin(x)? arent't we going in circles here?

• sin(a)/cos(a)/factorial(i) become known at compile-time so we can precalc them and hardcode them in the source
• the for-loop can be unrolled
• in fact, all the code "collapses" into "just another polynomial". that's because that's what this "taylor series thing" is: turning arbitrary functions into polynomials.

if you don't know what a polynomial is:

• it's a function of x like f(x) = c0 + c1*x + c2*(x**2) + c3*(x**3) + …
• c0,c1,c2,… are called coefficients. these are the constants that all this "taylor series stuff" helps us calculate
• the "…" means you can decide how many "parts" you want in the polynomial; each time you add a part, the "order of the polynomial" goes up by 1. f(x)=c0 is called a 0th order polynomial, f(x)=c0 + c1*x is a 1st order, etc.
• loosely speaking, the more orders a polynomial has, the more "approximating power" it has, while also being more expensive to calculate.

taylor series can also approximate other (but not all) non-siney functions, but for today we'll stick with sines. here's the general taylor series formula:

$$f\left(x\right)=f\left(a\right)+\frac{f\text{'}\left(a\right)}{1!}{(x-a)}^1+\frac{f\text{'}\text{'}\left(a\right)}{2!}{(x-a)}^2+\frac{f\text{'}\text{'}\text{'}\left(a\right)}{3!}{(x-a)}^3+\dots$$

• f(x) is the function we're trying to approximate
• a is the "anchor point", chosen by us
• f'(x) is the first derivative of f(x), f''(x) is the second, and so on. so f'(a) is the value of f'(x) for x=a.
• N! is the factorial of integer N:
  • 1! = 1
  • 2! = 2*1 = 2
  • 3! = 3*2*1 = 6
  • 4! = 4*3*2*1 = 24
  • (and you know what lets just define it recursively)
  • N! = N*(N-1)! for N>1 (use 1!=1 to terminate recursion)
• "..." means "continue the pattern if you like"; each time you add a section the order goes up by 1 (the order matches the order of the polynomial). a 0th order taylor series is just f(x)=f(a) (you can even try it in the plot above — open in new tab to avoid scrolling )

if you're new to derivatives they might seem a bit scary (sure scares me from time to time). the tool we're using is called differentiation, and it's a tool in the calculus-toolbox. when you "differentiate a function" f(x), you find its "derivative" f'(x). now you have stuff you can look up! but tl;dr: f'(x) is a function of x, just like f(x), except it returns the slope of f(x) (and f''(x) returns the slope of f'(x) and so on); try setting order=1 above; (open in new tab to avoid scrolling) this gives you a 1st order polynomial, also known as: "a line". notice how the line always tangents ("touches, but only barely") f(x) at x=a (drag the a-slider around to see it). that's because its slope is f'(a).

ok sure but how do you find the derivative of a function then? yeah well that's beyond the scope of this post, but i'll give you the cheatsheet for sin(x):

• f(x)=sin(x)
• f'(x)=cos(x)
• f''(x)=-sin(x)
• f'''(x)=-cos(x)
• f''''(x)=sin(x) (same as f(x))
• f'''''(x)=cos(x) (same as f'(x))
• and so on

to get a more intuitive feel for differentiation, since we've said that "f'(x) is the slope of f(x)", lets consider sin(x) in various places:

• around sin(0): it's the steepest that sin(x) gets (45°) so we expect the derivative to also be greatest there. sin'(0)=cos(0)=1, and cos(x) goes between [-1;1]. the math checks out. if we could only approximate sin(x) with a line, we could apply this slope: y=cos(0)*x or y=x. you can try this with order=1, a=0 in the plotter above, which gives you the y=x line.
• around sin(pi/2): the curve is flat here, so we expect the slope to be zero. cos(pi/2)=0. the math checks out. also, since sin(x) starts to drop for x>pi/2 we expect the derivative to drop below zero, and cos(x) crosses 0, into the negative, for x>pi/2. the math checks out.
• and so on. the math checks out.

you can find these and other differentiation cheatsheets on wikipedia and elsewhere. the syntax is different (of course) due to some tabs-vs-spaces kind of nerd war in the late 17th century, but d/dx sin(x) = cos(x) means the same as sin'(x) = cos(x) *ducks*

we can also briefly explore something called numerical differentiation, which is a generic way to differentiate a function without actually knowing any "differentiation rules".

to derive sin(x) we find the slope of the line between sin(x-delta) and sin(x+delta), for some small delta. and in the code above we do it recursively to find higher order derivatives. one nasty thing here is that if you lower the delta, you get better results early on, but it blows up faster.

if you've played around with the controls you may have noticed that sometimes adjacent orders look the same. for instance, if a=0, then the plots for order=5 and order=6 look the same. this is because some of the derivatives are zero (which also makes the polynomial coefficient zero) when a=0:

• f(0)=sin(0)=0
• f'(0)=cos(0)=1
• f''(0)=-sin(0)=0
• f'''(0)=-cos(0)=-1
• and so on

here's a hack to try in the editor above (open in new tab to avoid scrolling)

• set order=7
• set a=0
• set "reshape" checkmark
• verify that max error is around 0.00015
• in the bottom of bad(), change

  return sum;

  to

  return sum*1.00009;

if you need sin(pi/2) to be 1 as exactly as possible, you can also try "auto-calibrating" it:

• in the bottom of the code, outside of the bad() function, insert:

  print(1/bad(Math.PI/2))

• press ctrl+enter and copy the number from the debug output below your code (mine says 1.000156923221535)
• use that instead in the bad() return statement:

  return sum*1.000156923221535;

• now max error is about 0.00013, not much better, but we have the nice property that sin(0)=0, sin(pi/2)=1, etc

to get back to plotting, lets find a fine sine. to get the polynomial coefficients we can insert this in the bottom of the for-loop in bad():

print([i,fa/factorial(i)]);

bad(0);

so applying everything until now for an order=9 sine, with sin(0)=0, sin(pi/2)=0, i get...

in C this actually performs slightly worse than our previous best sine from last post (which puzzles me a little bit because the previous one has a division, and the same number of multiplications, but maybe it has better parallelism?). and in terms of accuracy it's also worse! document your failures! some say there are more to come.

but one cool thing you can't take away from it: we analyzed and plotted and understood what was going on (if only a little bit) instead of just copying coefficients from wikipedia.

we have a reusable tool now!

we're going to be using it soon to implement exp2(x). but before we get to that we need another tool: delicious low-level floating-point binary hackery!

also re: Horner's method: polynomials like

c0 + c1*(x**1) + c2*(x**2) + c3*(x**3) + ...

c0 + x*(c1 + x*(c2 + x*(c3 + ...)))

also-also a few last words on differentiation: it's the "magic part" of taylor series where you'll find its strengths and weaknesses. some software like wxmaxima and wolfram alpha can differentiate functions for you (try if it can turn sin(x) into cos(x)). in the code above, you only need to change the code that calculates fa in order to approximate a completely different function (and maybe remove the triangle wave stuff that's tailored to sines)

but taylor series can't deal with functions that have discontinuities, unless their derivatives somehow "know about them". consider a square wave, defined as (mixed math/JS syntax)

f(x)=(x%1)<.5?1:-1

f'(x)=0

however, froos has before talked about fourier series, which can approximate square waves and, well, anything, and these can be differentiated! probably still doesn't mean taylor series will be any good at sharp edges, but at least they'll have a chance.