Lesson 32: Iteration and approximation restart;

Recurrence relations are often used in numerical methods. You want to compute some sequence LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUYvNiVRIm5GJ0YyRjUvRjZRJ25vcm1hbEYnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPQ== which satisfies a recurrence relation, so you start with known values for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMEYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHRj1GPg== or the first few LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUYvNiVRIm5GJ0YyRjUvRjZRJ25vcm1hbEYnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPQ==, and iterate the recurrence formula. This may or may not be a good way to calculate LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUYvNiVRIm5GJ0YyRjUvRjZRJ25vcm1hbEYnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPQ==. The main thing that can go wrong is that the inevitable roundoff errors can grow with each iteration until they overwhelm the true solution. For example, suppose you want to approximate the following sequence of numbers, related to the remainders in the series for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiZUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Jy1JJW1zdWJHRiQ2JS1GLDYlUSJ3RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtRiw2JVEibkYnRjdGOi9GO1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RIj1GJ0ZCLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZNLyUpc3RyZXRjaHlHRk0vJSpzeW1tZXRyaWNHRk0vJShsYXJnZW9wR0ZNLyUubW92YWJsZWxpbWl0c0dGTS8lJ2FjY2VudEdGTS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRmZuLUYjNidGKy1GIzYlRj8tRkg2LVEiIUYnRkJGS0ZORlBGUkZURlZGWC9GZW5RLDAuMTExMTExMWVtRicvRmhuRmFvRkItRkg2LVExJkludmlzaWJsZVRpbWVzO0YnRkJGS0ZORlBGUkZURlZGWC9GZW5RJjAuMGVtRicvRmhuRmdvLUkobWZlbmNlZEdGJDYkLUYjNiYtRkg2LVEvJkV4cG9uZW50aWFsRTtGJ0ZCRktGTkZQRlJGVEZWRlhGZm9GYm8tRkg2LVEoJm1pbnVzO0YnRkJGS0ZORlBGUkZURlZGWC9GZW5RLDAuMjIyMjIyMmVtRicvRmhuRmVwLUZqbzYkLUYjNictSSttdW5kZXJvdmVyR0YkNictRkg2L1EmJlN1bTtGJy8lK2ZvcmVncm91bmRHUS5bMTQ0LDE0NCwxNDRdRidGQi9JK21zZW1hbnRpY3NHRiRRJmluZXJ0RidGS0ZOL0ZRRjlGUi9GVUY5L0ZXRjlGWEZmby9GaG5RLDAuMTY2NjY2N2VtRictRiM2Ji1GLDYlUSJrRidGN0Y6RkctSSNtbkdGJDYkRkZGQkZCRj9GWC8lLGFjY2VudHVuZGVyR0ZNRistSSdtc3BhY2VHRiQ2Ji8lJ2hlaWdodEdRJjAuMGV4RicvJSZ3aWR0aEdRJDUuMEYnLyUmZGVwdGhHRltzLyUqbGluZWJyZWFrR1ElYXV0b0YnLUkmbWZyYWNHRiQ2KC1GIzYkLUZicjYkUSIxRidGQkZCLUYjNiZGKy1GIzYlRl5yRl1vRkJGK0ZCLyUubGluZXRoaWNrbmVzc0dGW3QvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGZHQvJSliZXZlbGxlZEdGTUZCRkJGQkZCRkJGK0ZCRitGQg==. w:= n -> n! * (exp(1) - Sum(1/k!,k=0..n)); Digits:= 10: seq(evalf(w(n)),n=0..15); Catastrophic cancellation is giving us fewer and fewer significant digits, because we're subtracting two numbers LUkjbW9HNiMvSSttb2R1bGVuYW1lRzYiSSxUeXBlc2V0dGluZ0dJKF9zeXNsaWJHRic2LVEiZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGMS8lKXN0cmV0Y2h5R0YxLyUqc3ltbWV0cmljR0YxLyUobGFyZ2VvcEdGMS8lLm1vdmFibGVsaW1pdHNHRjEvJSdhY2NlbnRHRjEvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZA and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkrbXVuZGVyb3ZlckdGJDYnLUkjbW9HRiQ2L1EmJlN1bTtGJy8lK2ZvcmVncm91bmRHUS5bMTQ0LDE0NCwxNDRdRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnL0krbXNlbWFudGljc0dGJFEmaW5lcnRGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGPS8lKXN0cmV0Y2h5R1EldHJ1ZUYnLyUqc3ltbWV0cmljR0Y9LyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRj0vJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4xNjY2NjY3ZW1GJy1GIzYmLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR0ZCL0Y2USdpdGFsaWNGJy1GLzYtUSI9RidGNUY7Rj4vRkFGPUZDL0ZGRj0vRkhGPUZJL0ZMUSwwLjI3Nzc3NzhlbUYnL0ZPRlxvLUkjbW5HRiQ2JFEiMEYnRjVGNS1GVDYlUSJuRidGV0ZZRkkvJSxhY2NlbnR1bmRlckdGPS1GVDYjUSFGJy1JJ21zcGFjZUdGJDYmLyUnaGVpZ2h0R1EmMC4wZXhGJy8lJndpZHRoR1EkNS4wRicvJSZkZXB0aEdGX3AvJSpsaW5lYnJlYWtHUSVhdXRvRictSSZtZnJhY0dGJDYoLUYjNiQtRl9vNiRRIjFGJ0Y1RjUtRiM2JkZnby1GIzYlRlMtRi82LVEiIUYnRjVGO0Y+RmhuRkNGaW5Gam5GSS9GTFEsMC4xMTExMTExZW1GJy9GT0ZocUY1RmdvRjUvJS5saW5ldGhpY2tuZXNzR0ZfcS8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0Zeci8lKWJldmVsbGVkR0Y9RjU= that are very close. Can we get around that? We could increase Digits, but for any particular value of Digits we'll still have a problem when LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEibkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= is large enough (and arbitrary precision isn't always available). It's easy to see that we have the following recurrence relation: recurrence := y(n) = n*y(n-1) - 1; eval(recurrence, y=w); expand(%); simplify(lhs(%)-rhs(%)); combine(%); value(%); This suggests the following method. W[0]:= evalf(exp(1)-1); for count from 1 to 10 do W[count]:= count*W[count-1] - 1 end do; So far, so good? But... for count from 11 to 20 do W[count]:= count*W[count-1] - 1 end do; These results are absurd. The LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJW1zdWJHRiQ2JS1GLDYlUSJXRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtRiw2JVEibkYnRjdGOi9GO1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCRitGQg== should all be positive, and not very large. Roundoff error has struck again! What happened here is that whatever roundoff error is in LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJW1zdWJHRiQ2JS1GLDYlUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiZGKy1GIzYmLUYsNiVRIm5GJ0Y3RjotSSNtb0dGJDYtUSgmbWludXM7RicvRjtRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZMLyUpc3RyZXRjaHlHRkwvJSpzeW1tZXRyaWNHRkwvJShsYXJnZW9wR0ZMLyUubW92YWJsZWxpbWl0c0dGTC8lJ2FjY2VudEdGTC8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRmVuLUkjbW5HRiQ2JFEiMUYnRkhGSEYrRkgvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZIRitGSA== gets multiplied by n in LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJW1zdWJHRiQ2JS1GLDYlUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtRiw2JVEibkYnRjdGOi9GO1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCRitGQg==. So any error in LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJW1zdWJHRiQ2JS1GLDYlUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtSSNtbkdGJDYkUSIwRicvRjtRJ25vcm1hbEYnRkMvJS9zdWJzY3JpcHRzaGlmdEdGQkZDRitGQw== would be multiplied by LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiUtRiw2JVEibkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIiFGJy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4xMTExMTExZW1GJy8lJ3JzcGFjZUdGU0ZARitGQEYrRkA= by the time you get to LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJW1zdWJHRiQ2JS1GLDYlUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtRiw2JVEibkYnRjdGOi9GO1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCRitGQg==, even if no further roundoff error was introduced. Eventually the roundoff errors are larger than the actual values. How can we get around this? By using a different numerical method, one that doesn't magnify roundoff errors. In this case, one way is to use the recursion but go backwards: y(n-1) = solve(recurrence,y(n-1)); If we use an approximate value for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJuRidGNEY3Rj5GPkY+RitGPg== to calculate LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2JkYrLUYjNiYtRiw2JVEibkYnRjRGNy1GOzYtUSgmbWludXM7RidGPkZARkNGRUZHRklGS0ZNL0ZQUSwwLjIyMjIyMjJlbUYnL0ZTRlxvLUkjbW5HRiQ2JFEiMUYnRj5GPkYrRj5GPkY+RitGPg==, any error in the value for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJuRidGNEY3Rj5GPkY+RitGPg== gets divided by LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEibkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=. The errors decrease, so LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtRiw2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnL0Y6USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGRC8lKXN0cmV0Y2h5R0ZELyUqc3ltbWV0cmljR0ZELyUobGFyZ2VvcEdGRC8lLm1vdmFibGVsaW1pdHNHRkQvJSdhY2NlbnRHRkQvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZTLUkobWZlbmNlZEdGJDYkLUYjNiQtSSNtbkdGJDYkUSIwRidGQEZARkBGQEYrRkBGK0ZA should end up very accurate (of course, additional error is introduced by roundoff at each step, but it won't be magnified). We could start off with even a very rough approximation for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtRiw2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnL0Y6USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGRC8lKXN0cmV0Y2h5R0ZELyUqc3ltbWV0cmljR0ZELyUobGFyZ2VvcEdGRC8lLm1vdmFibGVsaW1pdHNHRkQvJSdhY2NlbnRHRkQvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZTLUkobWZlbmNlZEdGJDYkLUYjNiQtRiw2JVEiTkYnRjZGOUZARkBGQEYrRkBGK0ZA, and after a few backward steps the error will be quite small. W[50] := 0.0; for nn from 50 to 1 by -1 do W[nn-1] := (W[nn]+1)/nn end do; Actually LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Jy1JJW1zdWJHRiQ2JS1GLDYlUSJ3RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtRiw2JVEibkYnRjdGOi9GO1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RIj1GJ0ZCLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZNLyUpc3RyZXRjaHlHRk0vJSpzeW1tZXRyaWNHRk0vJShsYXJnZW9wR0ZNLyUubW92YWJsZWxpbWl0c0dGTS8lJ2FjY2VudEdGTS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRmZuLUYjNictSSttdW5kZXJvdmVyR0YkNictRkg2L1EmJlN1bTtGJy8lK2ZvcmVncm91bmRHUS5bMTQ0LDE0NCwxNDRdRidGQi9JK21zZW1hbnRpY3NHRiRRJmluZXJ0RidGS0ZOL0ZRRjlGUi9GVUY5L0ZXRjlGWC9GZW5RJjAuMGVtRicvRmhuUSwwLjE2NjY2NjdlbUYnLUYjNictRiw2JVEia0YnRjdGOkZHLUYjNiZGPy1GSDYtUSIrRidGQkZLRk5GUEZSRlRGVkZYL0ZlblEsMC4yMjIyMjIyZW1GJy9GaG5GaXAtSSNtbkdGJDYkUSIxRidGQkZCRitGQi1GLDYlUSgmaW5maW47RidGN0Y6RlgvJSxhY2NlbnR1bmRlckdGTUYrLUknbXNwYWNlR0YkNiYvJSdoZWlnaHRHUSYwLjBleEYnLyUmd2lkdGhHUSQ1LjBGJy8lJmRlcHRoR0ZpcS8lKmxpbmVicmVha0dRJWF1dG9GJy1JJm1mcmFjR0YkNigtRiM2JkYrLUYjNiVGPy1GSDYtUSIhRidGQkZLRk5GUEZSRlRGVkZYL0ZlblEsMC4xMTExMTExZW1GJy9GaG5GXXNGQkYrRkItRiM2JkYrLUYjNiVGYHBGaXJGQkYrRkIvJS5saW5ldGhpY2tuZXNzR0ZecS8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0Zncy8lKWJldmVsbGVkR0ZNRkJGK0ZCRitGQg==. Maple should be able to evaluate these rapidly converging series quite well. Let's compare its results to the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiV0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUYvNiVRIm5GJ0YyRjUvRjZRJ25vcm1hbEYnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPQ== we just computed. [seq([n,W[n]-evalf(Sum(n!/k!,k=n+1..infinity))],n=0..50)]; Up to n=45, our results were about as accurate as you could hope for with Digits = 10. If we wanted LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEid0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEjMjBGJy9GNlEnbm9ybWFsRidGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPg==, would it have been worthwhile to try to get a more accurate value for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEid0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEjNTBGJy9GNlEnbm9ybWFsRidGPi8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnRj4=? A method is numerically stable if small errors in intermediate results will have a negligible effect on the final result, or numerically unstable if they will have a large effect on the final result. So this recurrence was numerically unstable in the forward direction, but numerically stable in the backward direction. It's a good idea to look for numerical stability when choosing a numerical method.

Suppose you're interested in computing a particular function LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJmRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ4RidGNEY3Rj5GPkY+RitGPg==, perhaps a rather complicated one. You're going to need to compute it lots of times, so it's important to have an efficient way of doing that. For example, perhaps you're designing a device that will need to compute this function, and will need to do it very quickly. You're willing to spend some effort in finding an efficient way to find a good method of calculating the function, if that will be rewarded with improved performance when the device is operating. Fortunately, you don't need to compute it exactly: a good approximation will do. Specifically, you are given LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEtJnZhcmVwc2lsb247RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJ0Yy > 0 and an interval LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkobWZlbmNlZEdGJDYmLUYjNiYtSSNtaUdGJDYlUSJhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiLEYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GMTYlUSJiRidGNEY3Rj5GPi8lJW9wZW5HUSJbRicvJSZjbG9zZUdRIl1GJ0Y+, and you want an easily computed function LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJnRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ4RidGNEY3Rj5GPkY+RitGPg== such that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1JKG1mZW5jZWRHRiQ2KC1GIzYoRistRiM2Ji1GLDYlUSJnRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRj9RJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZJLyUpc3RyZXRjaHlHRkkvJSpzeW1tZXRyaWNHRkkvJShsYXJnZW9wR0ZJLyUubW92YWJsZWxpbWl0c0dGSS8lJ2FjY2VudEdGSS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlgtRjI2JC1GIzYkLUYsNiVRInhGJ0Y7Rj5GRUZFRkUtRkI2LVEoJm1pbnVzO0YnRkVGR0ZKRkxGTkZQRlJGVC9GV1EsMC4yMjIyMjIyZW1GJy9GWkZgby1GIzYmLUYsNiVRImZGJ0Y7Rj5GQUZlbkZFRitGRUZFL0krbXNlbWFudGljc0dGJFEkYWJzRicvJSVvcGVuR1EpJnZlcmJhcjtGJy8lJmNsb3NlR0ZccEZnby1GQjYtUSI8RidGRUZHRkpGTEZORlBGUkZUL0ZXUSwwLjI3Nzc3NzhlbUYnL0ZaRmNwLUYsNiVRLSZ2YXJlcHNpbG9uO0YnL0Y8RklGRUZFRitGRQ== for all LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= in this interval. For example, a polynomial LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJnRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ4RidGNEY3Rj5GPkY+RitGPg== with not too high a degree would be nice. The particular example I'll take is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2KEYrLUYjNiUtSSVtc3ViR0YkNiUtRiw2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMEYnL0Y9USdub3JtYWxGJ0ZFLyUvc3Vic2NyaXB0c2hpZnRHRkQtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ4RidGOUY8RkVGRUZFLUkjbW9HRiQ2LVEiPUYnRkUvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRlcvJSlzdHJldGNoeUdGVy8lKnN5bW1ldHJpY0dGVy8lKGxhcmdlb3BHRlcvJS5tb3ZhYmxlbGltaXRzR0ZXLyUnYWNjZW50R0ZXLyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGYG8tRiM2Ji1GLDYlUSdhcmN0YW5GJy9GOkZXRkUtRlI2LVEwJkFwcGx5RnVuY3Rpb247RidGRUZVRlhGWkZmbkZobkZqbkZcby9GX29RJjAuMGVtRicvRmJvRl1wLUZKNiQtRiM2JkYrLUYjNidGKy1GIzYmLUZCNiRRIjJGJ0ZFLUZSNi1RMSZJbnZpc2libGVUaW1lcztGJ0ZFRlVGWEZaRmZuRmhuRmpuRlxvRlxwRl5wRk5GRS1GUjYtUSIrRidGRUZVRlhGWkZmbkZobkZqbkZcby9GX29RLDAuMjIyMjIyMmVtRicvRmJvRmFxLUZCNiRRIjFGJ0ZFRkVGK0ZFRkVGRUYrRkVGK0ZF on the interval LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkobWZlbmNlZEdGJDYmLUYjNictSSNtaUdGJDYjUSFGJy1GIzYlLUkjbW9HRiQ2LVEqJnVtaW51czA7RicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y/LyUpc3RyZXRjaHlHRj8vJSpzeW1tZXRyaWNHRj8vJShsYXJnZW9wR0Y/LyUubW92YWJsZWxpbWl0c0dGPy8lJ2FjY2VudEdGPy8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRk4tSSNtbkdGJDYkUSIxRidGOkY6LUY3Ni1RIixGJ0Y6Rj0vRkFRJXRydWVGJ0ZCRkRGRkZIRkovRk1RJjAuMGVtRicvRlBRLDAuMzMzMzMzM2VtRidGUUY6RjovJSVvcGVuR1EiW0YnLyUmY2xvc2VHUSJdRidGOg== with LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUS0mdmFyZXBzaWxvbjtGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEiPUYnRjcvJSZmZW5jZUdGNi8lKnNlcGFyYXRvckdGNi8lKXN0cmV0Y2h5R0Y2LyUqc3ltbWV0cmljR0Y2LyUobGFyZ2VvcEdGNi8lLm1vdmFibGVsaW1pdHNHRjYvJSdhY2NlbnRHRjYvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZOLUklbXN1cEdGJDYlLUkjbW5HRiQ2JFEjMTBGJ0Y3LUYjNiUtRjs2LVEqJnVtaW51czA7RidGN0Y+RkBGQkZERkZGSEZKL0ZNUSwwLjIyMjIyMjJlbUYnL0ZQRmhuLUZVNiRRIjhGJ0Y3RjcvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnRjdGK0Y3. f0:= x -> arctan(2*x+1); plot(f0(x),x=-1 .. 1); Hearing the words "polynomial" and "approximation", a Calculus student might first think of Taylor series. However, these turn out not to be useful for our purpose. The idea of a Taylor series is that it gives a really good approximation of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUklbXN1YkdGJDYlLUkjbWlHRiQ2I1EhRictRiM2JEYuLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnRi4tRiM2JS1GLDYlLUYvNiVRImZGJy8lJ2l0YWxpY0dRJXRydWVGJy9GNVEnaXRhbGljRictRiM2JC1JI21uR0YkNiRGOUY0RjRGNy1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYvNiVRInhGJ0ZBRkRGNEY0RjRGLkY0 near one point. On the rest of the interval, however, the approximation might be terrible. taylor(f0(x), x=0, 11); TaylorApp:= unapply(convert(%, polynom), x); plot(f0(x) - TaylorApp(x), x=-1..1); The approximation is excellent for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= near 0, but poor near the ends of the interval, so the maximum error is about 2.8. In fact, the radius of convergence of this Taylor series is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkmbWZyYWNHRiQ2KC1GIzYkLUkmbXNxcnRHRiQ2Iy1JI21uR0YkNiRRIjJGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRidGNy1GIzYkRjNGNy8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGQS8lKWJldmVsbGVkR1EmZmFsc2VGJ0Y3, and for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkobWZlbmNlZEdGJDYoLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GNlEnbm9ybWFsRicvSSttc2VtYW50aWNzR0YkUSRhYnNGJy8lJW9wZW5HUSkmdmVyYmFyO0YnLyUmY2xvc2VHRj9GOkY4 larger than this the Taylor series is basically useless.

The next type of series to look at is a Fourier series. We're looking at a function on LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkobWZlbmNlZEdGJDYmLUYjNictSSNtb0dGJDYtUSomdW1pbnVzMDtGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjkvJSlzdHJldGNoeUdGOS8lKnN5bW1ldHJpY0dGOS8lKGxhcmdlb3BHRjkvJS5tb3ZhYmxlbGltaXRzR0Y5LyUnYWNjZW50R0Y5LyUnbHNwYWNlR1EsMC4yMjIyMjIyZW1GJy8lJ3JzcGFjZUdGSC1JI21uR0YkNiRRIjFGJ0Y0LUYxNi1RIixGJ0Y0RjcvRjtRJXRydWVGJ0Y8Rj5GQEZCRkQvRkdRJjAuMGVtRicvRkpRLDAuMzMzMzMzM2VtRidGS0Y0RjQvJSVvcGVuR1EiW0YnLyUmY2xvc2VHUSJdRidGNA==, so we need to relate the interval LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkobWZlbmNlZEdGJDYmLUYjNiktSSNtbkdGJDYkUSIwRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEiLEYnRjQvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHUSV0cnVlRicvJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRjg2LVEifkYnRjRGOy9GP0Y9RkFGQ0ZFRkdGSUZLL0ZPRk0tRjE2JFEiMkYnRjRGUS1JI21pR0YkNiVRJyYjOTYwO0YnLyUnaXRhbGljR0Y9RjRGNEY0LyUlb3BlbkdRIltGJy8lJmNsb3NlR1EiXUYnRjQ= that we used for Fourier series to LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkobWZlbmNlZEdGJDYmLUYjNictSSNtb0dGJDYtUSomdW1pbnVzMDtGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjkvJSlzdHJldGNoeUdGOS8lKnN5bW1ldHJpY0dGOS8lKGxhcmdlb3BHRjkvJS5tb3ZhYmxlbGltaXRzR0Y5LyUnYWNjZW50R0Y5LyUnbHNwYWNlR1EsMC4yMjIyMjIyZW1GJy8lJ3JzcGFjZUdGSC1JI21uR0YkNiRRIjFGJ0Y0LUYxNi1RIixGJ0Y0RjcvRjtRJXRydWVGJ0Y8Rj5GQEZCRkQvRkdRJjAuMGVtRicvRkpRLDAuMzMzMzMzM2VtRidGS0Y0RjQvJSVvcGVuR1EiW0YnLyUmY2xvc2VHUSJdRictRjE2LVEiLkYnRjRGN0Y6RjxGPkZARkJGREZUL0ZKRlVGNA== The appropriate saw function is saw := x -> x/Pi - 2*floor(x/(2*Pi)) - 1; plot(saw(x),x=-10 .. 10, discont=true); The periodic version of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMEYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHRj0tSShtZmVuY2VkR0YkNiQtRiM2JC1GLzYlUSJ4RidGMkY1Rj5GPkY+ is then LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JS1JJW1zdWJHRiQ2JS1GLDYlUSJmRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtSSNtbkdGJDYkUSIwRicvRjtRJ25vcm1hbEYnRkMvJS9zdWJzY3JpcHRzaGlmdEdGQi1JKG1mZW5jZWRHRiQ2JC1GIzYmRistRiM2Ji1GLDYlUSRzYXdGJ0Y3RjotSSNtb0dGJDYtUTAmQXBwbHlGdW5jdGlvbjtGJ0ZDLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZXLyUpc3RyZXRjaHlHRlcvJSpzeW1tZXRyaWNHRlcvJShsYXJnZW9wR0ZXLyUubW92YWJsZWxpbWl0c0dGVy8lJ2FjY2VudEdGVy8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRmBvLUZINiQtRiM2JC1GLDYlUSJ4RidGN0Y6RkNGQ0ZDRitGQ0ZDRkNGK0ZD. f1 := f0 @ saw; plot([f0(x/Pi-1), f1(x)],x=-10..10, colour=[red,blue]); Now for the Fourier series coefficients. Maple can't do these integrals symbolically. 1/Pi * int(f1(t)*cos(k*t),t=0..2*Pi) assuming k::integer; That LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEmZmxvb3JGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkobWZlbmNlZEdGJDYkLUYjNiYtRiw2I1EhRictRiM2Ji1JJm1mcmFjR0YkNigtSSNtbkdGJDYkUSIxRidGMi1GQzYkUSIyRidGMi8lLmxpbmV0aGlja25lc3NHRkUvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGTS8lKWJldmVsbGVkR0YxLUkjbW9HRiQ2LVExJkludmlzaWJsZVRpbWVzO0YnRjIvJSZmZW5jZUdGMS8lKnNlcGFyYXRvckdGMS8lKXN0cmV0Y2h5R0YxLyUqc3ltbWV0cmljR0YxLyUobGFyZ2VvcEdGMS8lLm1vdmFibGVsaW1pdHNHRjEvJSdhY2NlbnRHRjEvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0Zgby1GQDYoLUYjNiQtRiw2JVEidEYnL0YwUSV0cnVlRicvRjNRJ2l0YWxpY0YnRjItRiM2JC1GLDYlUSUmcGk7RidGL0YyRjJGSUZLRk5GUEYyRjpGMkYyLUZTNi5RIn5GJy8lMGZvbnRfc3R5bGVfbmFtZUdRJVRleHRGJ0YyRlZGWEZaRmZuRmhuRmpuRlxvRl5vRmFvRjI=should be 0 for t in this interval. It's curious that Maple can't simplify this: simplify(floor(t/(2*Pi))) assuming t>0,t<2*Pi; That should be 0. If we manually make it 0, will it work? Ak:=eval(%%, floor=0) assuming k::posint; simplify(%) assuming k::posint; What about for particular values of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=? eval(%,k=2); eval(%%,k=0); Except for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIjBGJ0Y5Rjk=, it still can't do the integral in closed form. No problem, we'll do the coefficients using numerical integration. Knowing this, it's better to use Int rather than int, to save us from wasting time with symbolic integrations that won't work. a:= unapply(1/Pi*Int(arctan(2*t/Pi-1)*cos(k*t),t=0..2*Pi), k); a(0):= value(a(0)); seq(evalf(a(k)),k=0..5); Similarly for b(k). b:= unapply(1/Pi * Int(arctan(2*t/Pi-1)*sin(k*t),t=0..2*Pi),k); For an approximation to the sum of the series, we can use a partial sum, i.e. take the sum up to some N instead of up to infinity. Psum := N -> evalf(a(0)/2 + add(a(k)*cos(k*t) + b(k)*sin(k*t),k=1..N)); Psum(5); plot([f1(t),Psum(5)],t=-10..10,colour=[red,blue]); Here's an animation to see how well the approximation goes as N increases. I'm using display(..., insequence=true) rather than animate to avoid premature evaluation problems. with(plots): display([seq(plot([f1(t),Psum(N)],t=-1..7,title=('N'=N)),N=1..20)],insequence=true); The approximation is rather bad near the jumps at LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbW5HRiQ2JFEiMEYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJ0Yv and 2LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEnJiM5NjA7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJ0Yy. That's inevitable: the partial sum is a continuous function, so it can't stay close to the discontinuous function LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPg== when that takes a jump. Somewhat surprising, perhaps, is the fact that the partial sum "overshoots" the jump, and that the size of the "overshoot", i.e. the difference between the maximum of Psum(N) and the maximum of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPg==, doesn't go to 0 as N increases. This is called the Gibbs phenomenon.

The poor approximation of our function by partial sums of its Fourier series is due to the discontinuities. In general, the smoother a function is, the faster its Fourier series coefficients go to 0, and so the better approximation the partial sums of the Fourier series are for the function. The best case is for an analytic function, where the coefficients go to 0 exponentially as LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEibkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSgmc3JhcnI7RidGOUY7Rj5GQEZCRkRGRkZIRkpGTUY1LUY2Ni5RKCZpbmZpbjtGJ0YvRjJGO0Y+RkBGQkZERkZGSEZKRk1GOQ==. So we might try replacing the saw function by something analytic that takes the interval LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkobWZlbmNlZEdGJDYmLUYjNiktSSNtbkdGJDYkUSIwRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEiLEYnRjQvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHUSV0cnVlRicvJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRjg2LVEifkYnRjRGOy9GP0Y9RkFGQ0ZFRkdGSUZLL0ZPRk0tRjE2JFEiMkYnRjRGUS1JI21pR0YkNiVRJyYjOTYwO0YnLyUnaXRhbGljR0Y9RjRGNEY0LyUlb3BlbkdRIltGJy8lJmNsb3NlR1EiXUYnRjQ= to LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkobWZlbmNlZEdGJDYmLUYjNictSSNtb0dGJDYtUSomdW1pbnVzMDtGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjkvJSlzdHJldGNoeUdGOS8lKnN5bW1ldHJpY0dGOS8lKGxhcmdlb3BHRjkvJS5tb3ZhYmxlbGltaXRzR0Y5LyUnYWNjZW50R0Y5LyUnbHNwYWNlR1EsMC4yMjIyMjIyZW1GJy8lJ3JzcGFjZUdGSC1JI21uR0YkNiRRIjFGJ0Y0LUYxNi1RIixGJ0Y0RjcvRjtRJXRydWVGJ0Y8Rj5GQEZCRkQvRkdRJjAuMGVtRicvRkpRLDAuMzMzMzMzM2VtRidGS0Y0RjQvJSVvcGVuR1EiW0YnLyUmY2xvc2VHUSJdRictRjE2LVEiLkYnRjRGN0Y6RjxGPkZARkJGREZUL0ZKRlVGNA== Well, the cosine function does that. plot(f0(cos(t)),t=-10..10); Now we have a nice analytic function, that will be very well approximated by its Fourier series. This has an additional, very important, benefit: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSRjb3NGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RidGNy8lJmZlbmNlR0Y2LyUqc2VwYXJhdG9yR0Y2LyUpc3RyZXRjaHlHRjYvJSpzeW1tZXRyaWNHRjYvJShsYXJnZW9wR0Y2LyUubW92YWJsZWxpbWl0c0dGNi8lJ2FjY2VudEdGNi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRk4tSShtZmVuY2VkR0YkNiQtRiM2JkYrLUYjNiYtRiw2JVEia0YnL0Y1USV0cnVlRicvRjhRJ2l0YWxpY0YnLUY7Ni1RMSZJbnZpc2libGVUaW1lcztGJ0Y3Rj5GQEZCRkRGRkZIRkpGTEZPLUYsNiVRInRGJ0ZlbkZnbkY3RitGN0Y3RjdGK0Y3 is a polynomial in LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSRjb3NGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RidGNy8lJmZlbmNlR0Y2LyUqc2VwYXJhdG9yR0Y2LyUpc3RyZXRjaHlHRjYvJSpzeW1tZXRyaWNHRjYvJShsYXJnZW9wR0Y2LyUubW92YWJsZWxpbWl0c0dGNi8lJ2FjY2VudEdGNi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRk4tSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ0RicvRjVRJXRydWVGJy9GOFEnaXRhbGljRidGN0Y3RjdGK0Y3, so our approximation becomes a polynomial in LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= rather than a trigonometric function. seq(expand(cos(k*t)), k = 2 .. 10); fp := unapply(f0(cos(t)),t); for count from 0 to 10 do a[count]:= evalf(1/Pi * Int(fp(t)*cos(count*t),t=-Pi..Pi)) end do: a[0]/2 + add(a[k]*cos(k*t), k = 1 .. 10); expand(%); eval(%,cos(t)=x); ChebApp:= unapply(%,x); The Cheb is short for "Chebyshev". This is called a Chebyshev series approximation for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JS1JJW1zdWJHRiQ2JS1GLDYlUSJmRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtSSNtbkdGJDYkUSIwRicvRjtRJ25vcm1hbEYnRkMvJS9zdWJzY3JpcHRzaGlmdEdGQi1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiVRInhGJ0Y3RjpGQ0ZDRkNGK0ZD on the interval [-1,1]. The polynomials LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1JJW1zdWJHRiQ2JS1GLDYlUSJURicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtRiw2JVEibkYnRjdGOi9GO1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnRkIvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRk0vJSlzdHJldGNoeUdGTS8lKnN5bW1ldHJpY0dGTS8lKGxhcmdlb3BHRk0vJS5tb3ZhYmxlbGltaXRzR0ZNLyUnYWNjZW50R0ZNLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGZm4tSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJzRidGN0Y6RkJGQkZCRitGQg== such that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2KEYrLUYjNiYtRiw2JVEkY29zRicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnRjkvJSZmZW5jZUdGOC8lKnNlcGFyYXRvckdGOC8lKXN0cmV0Y2h5R0Y4LyUqc3ltbWV0cmljR0Y4LyUobGFyZ2VvcEdGOC8lLm1vdmFibGVsaW1pdHNHRjgvJSdhY2NlbnRHRjgvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZQLUkobWZlbmNlZEdGJDYkLUYjNiZGKy1GIzYmLUYsNiVRIm5GJy9GN1EldHJ1ZUYnL0Y6USdpdGFsaWNGJy1GPTYtUTEmSW52aXNpYmxlVGltZXM7RidGOUZARkJGREZGRkhGSkZMRk5GUS1GLDYlUSJ0RidGZ25GaW5GOUYrRjlGOUY5LUY9Ni1RIj1GJ0Y5RkBGQkZERkZGSEZKRkwvRk9RLDAuMjc3Nzc3OGVtRicvRlJGZW8tRiM2Ji1JJW1zdWJHRiQ2JS1GLDYlUSJURidGZ25GaW4tRiM2JEZaRjkvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0Y8LUZUNiQtRiM2JkYrLUYjNiZGM0Y8LUZUNiQtRiM2JEZeb0Y5RjlGOUYrRjlGOUY5RitGOUYrRjk= are called the Chebyshev polynomials of the first kind. Let's see how well the Chebyshev series does at approximating LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMEYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHRj0tSShtZmVuY2VkR0YkNiQtRiM2JC1GLzYlUSJ4RidGMkY1Rj5GPi1JI21vR0YkNi1RIi5GJ0Y+LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZQLyUpc3RyZXRjaHlHRlAvJSpzeW1tZXRyaWNHRlAvJShsYXJnZW9wR0ZQLyUubW92YWJsZWxpbWl0c0dGUC8lJ2FjY2VudEdGUC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRmluRj4= plot(f0(x)-ChebApp(x),x=-1..1);

Actually Maple has commands to produce Chebyshev series approximations: chebpade and chebyshev in the numapprox package. chebpade stands for Chebyshev-Pad\303\251, and the Pad\303\251 part is because it can do something a bit more general than just a Chebyshev series. We can use chebpade to produce an approximation to a given degree N. with(numapprox): chebpade(f0(x),x=-1..1,10); The result is given in terms of T(n,x), which stand for the Chebyshev polynomials. Notice that the coefficient of each LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJURicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2Ji1GLDYlUSJuRidGNEY3LUY7Ni1RIixGJ0Y+RkAvRkRGNkZFRkdGSUZLRk1GTy9GU1EsMC4zMzMzMzMzZW1GJy1GLDYlUSJ4RidGNEY3Rj5GPkY+RitGPg== is the coefficient of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSRjb3NGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RidGNy8lJmZlbmNlR0Y2LyUqc2VwYXJhdG9yR0Y2LyUpc3RyZXRjaHlHRjYvJSpzeW1tZXRyaWNHRjYvJShsYXJnZW9wR0Y2LyUubW92YWJsZWxpbWl0c0dGNi8lJ2FjY2VudEdGNi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRk4tSShtZmVuY2VkR0YkNiQtRiM2JkYrLUYjNiYtRiw2JVEibkYnL0Y1USV0cnVlRicvRjhRJ2l0YWxpY0YnLUY7Ni1RMSZJbnZpc2libGVUaW1lcztGJ0Y3Rj5GQEZCRkRGRkZIRkpGTEZPLUYsNiVRInRGJ0ZlbkZnbkY3RitGN0Y3RjdGK0Y3 in our Fourier series. The actual definition of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJURicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2Ji1GLDYlUSJuRidGNEY3LUY7Ni1RIixGJ0Y+RkAvRkRGNkZFRkdGSUZLRk1GTy9GU1EsMC4zMzMzMzMzZW1GJy1GLDYlUSJ4RidGNEY3Rj5GPkY+RitGPg== as Chebyshev polynomials is contained in the orthopoly package. I'm not going to load that package, but instead use it this way: eval(%, T = orthopoly[T]); Notice that this is our ChebApp(x) (up to round-off error affecting the last digit of some coefficients). % - ChebApp(x); The maximum error in ChebApp(x) was about LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1JI21uR0YkNiRRIjRGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSNtb0dGJDYtUTEmSW52aXNpYmxlVGltZXM7RidGNS8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGPi8lKXN0cmV0Y2h5R0Y+LyUqc3ltbWV0cmljR0Y+LyUobGFyZ2VvcEdGPi8lLm1vdmFibGVsaW1pdHNHRj4vJSdhY2NlbnRHRj4vJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZNLUklbXN1cEdGJDYlLUYyNiRRIzEwRidGNS1GIzYlLUY5Ni1RKiZ1bWludXMwO0YnRjVGPEY/RkFGQ0ZFRkdGSS9GTFEsMC4yMjIyMjIyZW1GJy9GT0ZmbkYxRjUvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnRjVGK0Y1, still more than our LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUS0mdmFyZXBzaWxvbjtGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEiPUYnRjcvJSZmZW5jZUdGNi8lKnNlcGFyYXRvckdGNi8lKXN0cmV0Y2h5R0Y2LyUqc3ltbWV0cmljR0Y2LyUobGFyZ2VvcEdGNi8lLm1vdmFibGVsaW1pdHNHRjYvJSdhY2NlbnRHRjYvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZOLUklbXN1cEdGJDYlLUkjbW5HRiQ2JFEjMTBGJ0Y3LUYjNiUtRjs2LVEqJnVtaW51czA7RidGN0Y+RkBGQkZERkZGSEZKL0ZNUSwwLjIyMjIyMjJlbUYnL0ZQRmhuLUZVNiRRIjhGJ0Y3RjcvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnRjdGK0Y3. To attempt an approximation with error at most LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEtJnZhcmVwc2lsb247RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJ0Yy, we can try chebyshev. We give chebyshev our desired LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEtJnZhcmVwc2lsb247RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJ0Yy, and it will produce an approximation of whatever degree it needs. chebyshev(f0(x),x=-1..1,10^(-8)); The degree here will be 32.

numapprox package chebpade (in numapprox package) orthopoly[T] chebyshev (in numapprox package) LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=

LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=