Why use Fourier series instead of Taylor? In dynamical systems with linear differential equations, we almost always break up the function of independent variable in sines and cosines. 
But suppose that my function is smooth and periodic. Then what advantage do I get by using Fourier series instead Taylor?
Inside the radius of convergence, the Taylor series converges uniformly which is favourable, while the Fourier might only converge pointwise. 
So why use Fourier? I’d find using polynomials as more intuitive than sines or cosines.
 A: If you're asking about practical advantage, then you need to forget about everything you learned in analysis class. A practical calculation doesn't care whether something is pointwise or uniformly convergent. In fact, it is often very useful to use asymptotic series, which aren't even pointwise convergent. 
Broadly speaking, what makes a series useful is how numerically accurate a result you can get with it, while using only as many terms as is practical. "Rigorous" notions of convergence are not useful here because they talk about the limit of infinitely many terms, which obviously is never attained in practice. (For example, it is in principle true that $\cos(x)$ is described everywhere by its Taylor series, but try calculating $\cos(10^8)$ using that series and see how many terms you need to get a reasonable answer. The Taylor series is completely useless for this task.)
Fourier series are useful in this sense because many phenomena in nature exhibit spatial or temporal translational invariance. In the simplest cases, this renders problems diagonal in Fourier space, allowing you to write down the exact solution in one step. In more complicated cases, you can render the problem almost-diagonal in Fourier space, and treat it perturbatively. 
A: Many of the other answers are addressing the practicalities of expanding in Fourier series versus Taylor series. But there is at least one physical reason for choosing one over the other, and that is that the expansion coefficients of a vector written in an orthonormal basis reveal particular types of physical information about the system being described by the function, and the type of physical information that is revealed depends on the choice of basis:
The physical relevance of expansions in an orthonormal basis
This is directly related to some of the answers here that cast Fourier series in the language of linear algebra, where the sines and cosines compose an orthonormal basis for the space of functions you are looking at. In this context, an orthonormal basis has four nice properties.  The first is merely that orthonormal bases are nice calculationally in that they allow you to easily compute expansion coefficients by taking the inner product of the vector with a basis element. The other three are more physically relevant:


*

*The expansion coefficients in an orthonormal basis have a direct physical interpretation. In the context of basic mechanics, for instance, the $v_x$ attached to $\hat{x}$ in a velocity vector really is the velocity in the $x$ direction, in the following sense. Imagine moving along the $x$-axis at velocity $v_x$ (of the object) and watching the object as it moves: the object will recede from you but keep pace with you along that direction. So, importantly, the components of a vector in an orthonormal basis reveal physically relevant information about the system described by that vector.

*The type of physical information that's revealed depends on the orthonormal basis that's chosen. If you choose to represent a velocity vector in the basis $\{\hat{x},\hat{y},\hat{z}\}$, then you get the velocities in these directions (in the sense explained in the previous bullet point). If instead you choose the basis where one of the basis vectors lies along $\vec{v}$, then the component along this basis vector is the speed of the particle. Different bases reveal different information and conceal different information.

*Finally, the expansion coefficients in an orthonormal basis are unique, which is what really allows us to interpret the quantities physically. (What would the $x$-component of velocity really mean if it could take on two different values?)
The physical interpretation of expansion coefficients in a Fourier series
In the context of Fourier series, these bullet points have the following meaning. 


*

*The expansion coefficients tell us how much the function looks like a particular sine or cosine function. Let's consider one of the terms in the expansion, given by
$$a_n\cos\left(\omega_n t\right) + b_n\sin\left(\omega_n t\right),$$
where $\omega_n=2\pi n/T$. If we re-write this as
$$A_n\sin\left(\omega_n t+\phi_n\right),$$
where $A_n=\sqrt{a_n^2+b_n^2}$ and $\phi_n=\tan^{-1}(a_n/b_n)$, then we can interpret the coefficients as giving us two pieces of information. The amplitude $A_n$ tells us how much this particular function comes into the expansion, and therefore how much the original function looks like this particular sine function. In other words, $A_n$ tells us to what extent the original function oscillates with frequency $n/T$. In other words, the expansion coefficients are giving us the frequency content of the original signal, where the expansion coefficients quantify this relationship. (The phases $\phi_n$ would be relevant in, for instance, interference experiments, but for the purposes of this discussion, it's the amplitude that really matters.)

*A Fourier series (and a Fourier transform) yields the frequency (or wavelength, depending on the context) content of the function that is describing some physical quantity in some physical system. One can use other sets of orthogonal functions, too. For instance, in electrostatics, we can expand the electrostatic potential using Legendre polynomials, and these functions yield the$^1$ multipole (monopole, dipole, quadrupole, etc.) structure of the electrostatic potential. Different orthonormal bases reveal different types of physical information about the system.

*Uniqueness again allows us to do the above. If the Fourier series was actually a Fourier series$^2$, then this wouldn't really work.
Taylor series
As far as I understand, there is no "natural" inner product on the set of polynomials that makes the set of monomials $x^n$ orthonormal, and so we can't really interpret a Taylor series as an expansion in an orthonormal basis. However, we can$^3$ interpret a Taylor series as an expansion in a non-orthogonal basis, and so we get back a little of the interpretational power.
In particular, the expansion coefficients reveal structural information about the function, in that it reveals how constant the function is, how linear the function is, how quadratic, how cubic, and so on. This can have important physical consequences as well, as evidenced by the fact that the simple harmonic oscillator is ubiquitous: if a potential energy is "very quadratic", then the behavior of the system is very much like a simple harmonic oscillator.
This is not the type of information that is revealed via a Fourier series$^4$, and hence different bases reveal different information.
Numerical approximation
Finally, as mentioned in other answers, one or the other can aid numerical computations, but I would argue that this also has physical implications (or is at least informed by the physics). As an archetypical example, just see the previous section about Taylor series. In some sense, the approximation of the system as being a simple harmonic oscillator is a numerical approximation, and it is one that has both physical meaning and physical justification.

1. All but the first non-zero term in the expansion actually depend on the origin chosen, so using "the" here is not correct, and there's a little sophistry here in the sense that I'm stressing the importance of the uniqueness of the expansion. The point is still a good one, though.
2. We can of course choose different intervals over which to compute a Fourier series for a particular function, and in that sense a Fourier series is not unique. However, the interval is usually determined by the physical system we are describing, and so in that sense we get uniqueness back. In fact, sometimes we are actually computing a Fourier transform, in which case the interval is the entire real line.
3. In physics we can, anyway, to get an intuitive feel for what we're doing. And perhaps we can in mathematics, too, to the extent that they allow us. I don't know to what extent mathematicians think of Taylor series in this way, though.
4. Except that, as mentioned in another answer, there is a relationship between Fourier series and Taylor series that is revealed in the context of complex analysis. So there's another bit of sophistry in the name of getting my point across.
A: One of the deepest reasons to use the Fourier expansion is the idea of Harmonic Analysis induced by symmetries. The main topic is called Representation Theory and it forms the foundation to large parts of physics. 
Here's the idea: consider a (finite - or compact, for which many similar properties hold) group $G$ and a vector space $V$ on which it operates. Assume for simplicity the field of the vector space has characteristic 0 or at least not 2 (the latter case is the nastiest, with many exceptions). If confused, just assume we have an $\mathbb{R}$- or $\mathbb{C}$-based vector space.
This is to mean that group elements $g\in G$ act as vector space homomorphisms on the vector space $V$. I.e. any group element $g:V\to V$ can be considered a linear map over $V$ whose group multiplication respects homomorphism composition, i.e. $(g_1 \cdot g_2)(v) = (g_1 \circ g_2) (v)$ where $\cdot$ denotes the multiplications as group members while $\circ$ denotes the composition as homomorphisms.
The ensuing structure is called a group representation or $G$-space which has a lot of nice properties. Amongst these is the fact that, for a given group $G$, the $G$-space can be decomposed in only very specific ways into $G$-subspaces. They decompose uniquely into so-called irreducible $G$-subspaces which cannot be further decomposed. Of these irreducible subspaces, there exist only very specific ones for the given group $G$ (for finite groups, there are, in fact, only a finite number of inequivalent irreducible $G$-space types). 
It means that, no matter how your $G$-space looks, the decomposition will always contain the same irreducible types, albeit with different contributions - think of this as a kind of "prime-factor decomposition" of your $G$-space; its components will only depend on $G$. You can establish that before computing anything else. This is what molecular dynamics scientists used to establish degeneracy of excitations, purely based on symmetry considerations - before even computing the precise frequencies or knowing the dynamical parameters of the system.
If you have additional properties of your $V$, such as a scalar product, and the group operation respects that, you get additional nice properties, such as that elements from distinct irreducible $G$-subspaces of your $G$-space will be orthogonal (Schur's Lemma). In addition, if your group is abelian, i.e. commutative, and your field is $\mathbb{C}$, the irreducible spaces have dimension 1.
And now we come to the Fourier analysis: the Fourier components are nothing else than decomposition of your vector into its irreducible components when your group $G$ are the translations.
For simplicity, consider the group of translation operators $T_k$ (for technical reasons, assume a finite one-dimensional chain and periodic boundary conditions) which transforms $j \in [n]$ to $T_k(j) := (j+k-1) \mod n + 1$ (i.e. shift by $k$ and wrap around to 1 if past $n$). Given a (complex-valued) vector $v = (v_1, v_2, \dots, v_n)$, a translation operator $T_k$ will transform $v$ as per $T_k v := (v_{T_k^{-1}1}, v_{T_k^{-1}2}, \dots, v_{T_k^{-1}n})$, i.e. all coordinates of $v$ are shifted by $k$ and wrapped around. According to Representation Theory, this $G$-space with the $T_k$ forming the group can be decomposed into one-dimensional irreducible subspaces. It turns out that these are precisely the Fourier components, as is easily verifiable.
Namely: consider our problem's $G$-subspace $F_m := \langle (e^{1 \cdot im}, e^{2\cdot im},\dots,e^{n\cdot im}) \rangle$ (the angles indicate the space generated by the vector inside). The $i$ is the imaginary unit here, and $m$ selects the irreducible subspace (in your usual lingo, the harmonic). Clearly, this is irreducible, as it is one-dimensional. It is also indeed a $G$-space. It is enough to apply $T_k$ to $(e^{1 \cdot im}, e^{2\cdot im},\dots,e^{n\cdot im})$. This gives you $(e^{(1-k)im}, e^{(2-k)im},\dots,e^{(n-k)im})$, but that's the same as $e^{-ikm} \cdot (e^{1 \cdot im}, e^{2\cdot im},\dots,e^{n\cdot im})$, thus just a constant times the original vector and therefore in $F_m$. We have shown that $F_m$ is an irreducible $G$-space for the translations.
Since we have a scalar product, we automatically know from representation theory that the $F_m$ are orthogonal to each other for different $m$. This holds also for other relevant groups. So, for instance if you have the rotation groups on the sphere, you get the spherical harmonics as irreducible spaces. Again, no need to run through tedious integrations on the surface of the sphere - it pops out for all of these, automatically, through Schur's Lemma.
TL;DR: If you have a group which operates on a vector space, with a few additional conditions (group is finite or compact, though some results work for non-compact groups, too; field is nicely behaved, no character 2, for instance), you get a large set of guarantees as to how your vector space can be decomposed and how this decomposition will look like. Additional structure on the vector space, such as a scalar product etc., will lead to additional guarantees. The Fourier decomposition is the decomposition arising with the special case of the group being the translations of the coordinates and can easily be generalised. The relevant field is called "Representation Theory" and is one of the most elegant topics of mathematics with a strong influence on contemporary physics.
A: *

*The complex exponentials are eigenfunctions of the derivative and integral operators. So if you're analyzing linear differential equations, and using Fourier series, then you can consider each term on its own. If you use Taylor series you have to consider interactions between one term and other terms in the series. (This is also why we often write our Fourier series in terms of complex exponentials rather than sines and cosines)

*Extrapolation. If I have a function $f(x)$ and I approximate it on a region $[x_1, x_2]$ with a finite-length Taylor series $F_T(x)$, then outside of $[x_1, x_2]$ the Taylor series will tend to go to infinity. If I approximate it with a finite-length Fourier series, the series will remain bounded as $x\to\infty$.
A: Great question. One reason that complex exponential expansions (which end up turning on sines and cosines for real-valued problems) are more natural that Taylor series expansions is that they don't require picking a special point to expand around. In many situations the differential equation is translationally invariant, and there's no natural point to Taylor expand around, so you need to pick an arbitrary point. A general pattern in physics is that if your problem setup has some symmetry, you definitely want to take advantage of that symmetry in solving it.
Another issue is that, as The Photon mentioned, polynomials inevitably get unboundedly large at large $x$, which doesn't match up with the periodic nature of the solutions. For any finite-order Taylor expansion, you need to manually truncate the solution outside of a single fundamental period, which is a little awkward.
But probably the most important reason is that you are dealing with differential equations, and sines and cosines have the very special property of remaining unchanged (up to a scaling factor) after two derivatives (and the complex exponential version is unchanged up to a scaling factor after even a single derivative). As you gain experience with Fourier transforms, you'll see that this fact allows you to convert many linear differential equations into algebraic ones that are much easier to deal with. By contrast, diffentiating a polynomial takes you down the ladder to a lower-order polynomial, so you never get back to where you started, no many how many derivatives you take. This fact prevents you from taking advantage of that technique.
A: From the perspective of complex numbers, Fourier series are polynomials — or at least, Taylor series. For example, take the real period $1$ function $f$ with $f(x) = x$ for $-1/2 \le x < 1/2$. The Fourier series of $f$ is
$$f(x) = \sum_{n=1}^\infty \frac{(-1)^{n-1}}{\pi n} \sin 2\pi n x.$$
Using complex numbers, this can be proved as follows: for $-1/2 < x < 1/2$ we have
$$ f(x) = \Im \sum_{n=1}^\infty \frac{(-1)^{n-1}}{\pi n} \mathrm{e}^{2\pi \mathrm{i} n x} = \Im \frac{1}{\pi} \log(1 + \mathrm{e}^{2 \pi \mathrm{i} x}) = \Im \frac{1}{\pi} \mathrm{i} (\pi x) = x $$
where the final step uses that if $\phi$ is the argument of $1+\mathrm{e}^{2 \pi \mathrm{i} x}$ then $\tan \phi = \frac{\sin 2\pi  x}{1+ \cos 2 \pi x} = \tan \pi x$, and so $\phi = \pi x$ by our assumption on $x$. Since the Fourier series for $f$ is clearly periodic with period $1$, we have equality for all $x\not\in \mathbb{Z} + \frac{1}{2}$.
To make the connection with Taylor series, consider 
$$F(z) = \frac{\log(1+z)}{\pi} = \sum_{n=1}^\infty \frac{(-1)^{n-1}}{\pi n} z^n.$$
We have seen that $F(\mathrm{e}^{2\pi i x})$ has imaginary part $f(x)$. So $f$ is the imaginary part of the restriction of $F$ to the boundary of the unit circle, and the Fourier series of $f$ is simply the Taylor series of $F$, evaluated on the unit circle. This Taylor series converges for all $z$ with $|z| \le 1$ and $z \not= -1$, and correspondingly, the Fourier series converges for all $x \not\in \mathbb{Z} + \frac{1}{2}$.
A: Fourier series let you model your space of potential solutions as an infinite dimensional analogue of a finite dimensional Euclidean vector space. All of your formal analytical arguments can take place in that infinite dimensional space, where you can use your Euclidean intuition.
For example, in place of an inner product $\langle u,v\rangle = \sum_{i=1}^n u_i \cdot v_i$ of Euclidean vectors, you have the inner product $\langle f,g \rangle = \frac{1}{2\pi} \int_0^{1} f(t) g(t) \, dt$ of functions. Furthermore, in place of the "standard orthonormal basis" $e_1,...,e_n$ of Euclidean $n$-space, you have the "standard orthonormal bases" $e^{2\pi n t}$, $n \in \mathbb{Z}$, and this leads to Plancherel's Formula for the inner product. Of course, one still has to do the analysis to prove convergence of the integrals and sums which occur in place of the finite sums of Euclidean geometry.
This has many very practical applications. For example, the methods for solving linear differential equations can be understood (and in fact developed) using tools of linear algebra, in analogy to linear algebra carried out in finite dimensional Euclidean vector spaces. 
