Polynomial order of an approximation of a section of sine and numerical accuracy

I was playing with the idea that a sine function is periodic. But even within one period there are symmetries, namely the second fourth of a period is the mirror image of the first fourth and the second half is minus the first half. So if I could approximate the sine function from $0$ to $\pi/2$, then I would have enough to approximate it everywhere for real values.

For the approximation I wanted to use a polynomial of order $N$, and I wanted find it such that the polynomial minimizes the error between it and sine between $0$ and $\pi/2$. Therefore I came up with the following problem formulation:

$$\min_{\vec{a}} \int_0^{\frac{\pi}{2}} {\left(\sin(x) - \sum_{n=0}^N{a_n x^n}\right)^2 dx}. \tag{1}$$

This can be solved mostly analytically by setting the derivative of $(1)$ with respect to each $a_n$ equal to zero:

$$\int_0^{\frac{\pi}{2}} {x^m \sum_{n=0}^N{a_n x^n}dx} = \int_0^{\frac{\pi}{2}} {x^m \sin(x) dx}, \quad \forall\; m=0,1,\dots, N \tag{2}$$

$$\sum_{n=0}^N{\frac{a_n \left(\frac{\pi}{2}\right)^{n+m+1}}{n+m+1}} = \left\{\begin{array}{l,l} m!\left(\frac{\pi}{2}\right)^{m-1}\sum_{n=0}^{\frac{m-1}{2}}\frac{\left(\frac{-4}{\pi^2}\right)^n}{(m-2n-1)!}, & \text{if } m \text{ odd} \\ m!\left((-1)^{\frac{m-2}{2}}+\left(\frac{\pi}{2}\right)^{m-1}\sum_{n=0}^{\frac{m-2}{2}}\frac{\left(\frac{-4}{\pi^2}\right)^n}{(m-2n-1)!}\right), & \text{else}\end{array}\right.\tag{3}$$

Solving this problem then turns into solving a matrix equation ($Ax=b$).

But here is my question; if I increase $N$ I would expect getting smaller and smaller errors up to some point when the computational errors of calculating the polynomial coefficients get to big, mainly probably due to the sum of the right hand side of $(3)$. But how could these computational errors be decreased? Namely initially the error get smaller, but beyond a polynomial of order 8 the error seems to grow exponentially as can be seen in the following figure.

The error of the 8th order polynomial is at most $2.152\times 10^{-8}$ (while an 7th/8th order Taylor polynomial at $x=0$ has maximum error of $1.569\times 10^{-4}$). I just wonder if I could do something to increase the order and decrease the error to something closer to $10^{-15}$, but even using the variable-precision arithmetic in MATLAB does not seem to affect the error at all.

• Even for small $N=9$ the lhs matrix $A$ has condition number $3.8\times 10^{14}$, making this problem extremely ill-conditioned, which also explains why the error starts increasing with $N$. It is probably easiest to just pick a better basis for the least squares than the monomial basis, which is very poorly behaved. – Kirill May 8 '16 at 0:35
• @Kirill Might $\{x^N,x^{N-1}(x-\pi/2),\dots,x(x-\pi/2)^{N-1},(x-\pi/2)^N\}$ be a better choice for a basis? – fibonatic May 8 '16 at 1:02
• @fibonatic Informally, the shapes of the basis functions are too similar, which makes the optimal coefficients susceptible to small perturbations in the input data (having the same shape would make the coefficients completely ill-defined). I haven't checked, but that doesn't look like a good basis for the same reason. I agree with the Chebyshev interpolation suggestion, it is very well-behaved, especially on such simple problems; it would also make sense to consider some other loss function than the $L^2$-norm. – Kirill May 8 '16 at 1:43
• @Kirill What about orthonormal basis, such that the coefficients can simply be found by calculating the inner product (or does this not minimize the maximum error?). I am not familiar with interpolation at Chebyshev nodes, but from what I have read and tried so far they do not seem to be orthonormal. – fibonatic May 8 '16 at 12:30
• @fibonatic If you mean Legendre polynomials, those should work just fine. In fact, they are behaved very similarly to Chebyshev polynomials (they are all asymptotically similar to sines/cosines). Chebyshev are also orthogonal, but under a different choice of weight function in the norm. The reason for suggesting Chebyshev is that if you replace the optimization objective with the condition that weighted inner product with a Chebyshev polynomials is zero, that reduces the problem to a DCT, which is quite straightforward. But any orthogonal system will be much better conditioned. – Kirill May 8 '16 at 16:13

By using orthonormal polynomial basis functions, then you do not need to solve a linear set of equations in order to get the polynomial coefficients, similar to the Fourier series coefficients. Because for large $N$ the matrix, of the linear set of equations, will get ill-conditioned.

It is at least required that the basis functions are orthogonal, but normalized functions help with the ease of approximation of a function. In order to get slightly simpler basis functions (without transcendental coefficients induced by the upperbound of the interval of $\pi/2$) the interval will be changed into $[0,1]$, such that the function we want to approximate becomes $\sin\left(\frac{\pi}{2}x\right)$.

Two functions are orthogonal when there inner product is equal to zero. This inner product for this interval can be written as,

$$\left<f(x), g(x)\right> = \int_0^1 {f(x) g(x) dt}, \tag{1}$$

and a functions $f(x)$ is normalized when $\left<f(x), f(x)\right>=1$.

In order to get similar converges of the approximation each next orthonormal polynomial basis function will contain a higher power of $x$. The first normalized basis function of order zero will then become,

$$p_0(x) = 1. \tag{2}$$

Each following orthonormal polynomial basis function can then be found with,

$$p_k^*(x) = x^k - \sum_{i=0}^{k-1}\left<x^k, p_i(x)\right> p_i(x), \tag{3a}$$

$$p_k(x) = \frac{p_k^*(x)}{\sqrt{\left<p_k^*(x), p_k^*(x)\right>}}. \tag{3b}$$

The next three basis functions for example then become,

$$\begin{array}{} p_1(x) = \sqrt{3} \left(1 - 2 x\right) \\ p_2(x) = \sqrt{5} \left(1 - 6 x + 6 x^2\right) \\ p_3(x) = \sqrt{7} \left(1 - 12 x + 30 x^2 - 20 x^3\right) \end{array} \tag{4}$$

The $n$th order approximation of a function $f(x)$ on the interval $[0,1]$ can then simply be found with,

$$f_n(x) = \sum_{i=0}^n {\left<f(x), p_i(x)\right> p_i(x)}. \tag{5}$$

In the following figure it can be seen that the errors keeps decreasing for higher order of approximation of $\sin\left(\frac{\pi}{2}x\right)$.

• Aren't the $p_k(x)$ just the Legendre polynomials $P_k(2x-1)$ times a normalizing factor of $\sqrt{2k+1}$? $x\mapsto 2x-1$ maps the interval $[0,1]$ to $[-1,1]$, so to get orthonormal polynomials on $[0,1]$ or $[0,\frac\pi2]$ it's enough to map the interval to $[-1,1]$ by a linear change of variable. [BTW, please accept your own answer, stackoverflow.com/help/self-answer ] – Kirill May 10 '16 at 15:35
• @Kirill Your are correct that these are indeed remapped Legendre polynomials. The only down side of these are that for high order the absolute functional value at the boundary of the interval becomes large (indeed $\sqrt{2k+1}$), so the convergence at the edges of the boundary is not that great. Maybe this can be reduced by choosing a bigger interval ($\sin\left(\frac{\pi}{2}((1+2\beta)x-\beta\right)$). Also by looking at the error plot of my question it seems that the solutions are the same (before the matrix becomes too ill-conditioned). Is there an easy way to proof this? – fibonatic May 10 '16 at 16:24
• The polynomial $p(x)$ of a given degree that minimizes $\|f(x) - p(x)\|^2$ is unique, so whether you express it through its coefficients in the monomial basis or the Legendre basis, it still has the same values. The coefficients should be $G^{-1}\langle f, \phi_m\rangle$ where $G$ is the Gram matrix of the basis $\{\phi_m\}_{m}$ (it's diagonal for orthogonal polynomials). – Kirill May 10 '16 at 16:46