| 4 | |
| 5 | == Minimax polynomial for sin(x) == |
| 6 | |
| 7 | === Absolute error === |
| 8 | |
| 9 | Suppose we want to approximate sin(x) on [-π/2; π/2] with a polynomial P(x) such that the '''absolute''' error is never more than E: |
| 10 | |
| 11 | {{{ |
| 12 | #!latex |
| 13 | \[\big\vert\sin(x) - P(x)\big\vert \le E \qquad \forall x \in \bigg[-\frac{\pi}{2}, \frac{\pi}{2}\bigg]\] |
| 14 | }}} |
| 15 | |
| 16 | We know sin(x) is an odd function, so instead we look for a polynomial Q(x) such that P(x) = xQ(x²): |
| 17 | |
| 18 | {{{ |
| 19 | #!latex |
| 20 | \[\big\vert\sin(x) - xQ(x^2)\big\vert \le E \qquad \forall x \in \bigg[-\frac{\pi}{2}, \frac{\pi}{2}\bigg]\] |
| 21 | }}} |
| 22 | |
| 23 | Substitute y for x² and reduce the range to positive values: |
| 24 | |
| 25 | {{{ |
| 26 | #!latex |
| 27 | \[\big\lvert\sin(\sqrt{y}) - \sqrt{y}Q(y)\big\rvert \le E \qquad \forall y \in \bigg[0, \frac{\pi^2}{4}\bigg]\] |
| 28 | }}} |
| 29 | |
| 30 | Divide through by √y: |
| 31 | |
| 32 | {{{ |
| 33 | #!latex |
| 34 | \[\bigg\lvert\frac{\sin(\sqrt{y})}{\sqrt{y}} - Q(y)\bigg\rvert \le \frac{E}{|\sqrt{y}|} \qquad \forall y \in \bigg[0, \frac{\pi^2}{4}\bigg]\] |
| 35 | }}} |
| 36 | |
| 37 | If we want to force the asymptotic behaviour in x=0, we substitute Q(y) with 1+yR(y): |
| 38 | |
| 39 | {{{ |
| 40 | #!latex |
| 41 | \[\bigg\lvert\frac{\sin(\sqrt{y})}{\sqrt{y}} - 1 - yR(y)\bigg\rvert \le \frac{E}{|\sqrt{y}|} \qquad \forall y \in \bigg[0, \frac{\pi^2}{4}\bigg]\] |
| 42 | }}} |
| 43 | |
| 44 | Divide through by y: |
| 45 | |
| 46 | {{{ |
| 47 | #!latex |
| 48 | \[\bigg\lvert\frac{\sin(\sqrt{y})-\sqrt{y}}{y\sqrt{y}} - R(y)\bigg\rvert \le \frac{E}{|y\sqrt{y}|} \qquad \forall y \in \bigg[0, \frac{\pi^2}{4}\bigg]\] |
| 49 | }}} |
| 50 | |
| 51 | We then use the following code: |
| 52 | |
| 53 | {{{ |
| 54 | #!cpp |
| 55 | |
| 56 | static real myfun(real const &x) |
| 57 | { |
| 58 | real y = sqrt(x); |
| 59 | return (sin(y) - y) / (x * y); |
| 60 | } |
| 61 | |
| 62 | static real myerr(real const &x) |
| 63 | { |
| 64 | real y = sqrt(x); |
| 65 | return re(x * y); |
| 66 | } |
| 67 | |
| 68 | RemezSolver<6> solver; |
| 69 | solver.Run(real::R_1 >> 400, real::R_PI_2 * real::R_PI_2, myfun, myerr, 15); |
| 70 | }}} |
| 71 | |
| 72 | These are the resulting R, Q and P: |
| 73 | |
| 74 | {{{ |
| 75 | #!latex |
| 76 | \begin{eqnarray*} |
| 77 | R(y) & = & a_0 + a_1 y + a_2 y^2 + a_3 y^3 + a_4 y^4 + a_5 y^5 + a_6 y^6 \\ |
| 78 | Q(y) & = & 1 + a_0 y + a_1 y^2 + a_2 y^3 + a_3 y^4 + a_4 y^5 + a_5 y^6 + a_6 y^7 \\ |
| 79 | P(x) & = & x + a_0 x^3 + a_1 x^5 + a_2 x^7 + a_3 x^9 + a_4 x^{11} + a_5 x^{13} + a_6 x^{15} \\ |
| 80 | a_0 & = & -1.66666666666658080941942898789420724.10^{-1} \\ |
| 81 | a_1 & = & +8.33333333326271609442503773834687308.10^{-3} \\ |
| 82 | a_2 & = & -1.98412698200591143928364634696492885.10^{-4} \\ |
| 83 | a_3 & = & +2.75573160733868922065738227278330896.10^{-6} \\ |
| 84 | a_4 & = & -2.50518513021429359590028300127165228.10^{-8} \\ |
| 85 | a_5 & = & +1.60472959182597740337401201006549498.10^{-10} \\ |
| 86 | a_6 & = & -7.36458957326227991327065122848667046.10^{-13} \\ |
| 87 | \end{eqnarray*} |
| 88 | }}} |