Why do we care about $L^p$ spaces besides $p = 1$, $p = 2$, and $p = \infty$? I was helping a student study for a functional analysis exam and the question came up as to when, in practice, one needs to consider the Banach space $L^p$ for some value of $p$ other than the obvious ones of $p=1$, $p=2$, and $p=\infty$.  I don't know much analysis and the best thing I could think of was Littlewood's 4/3 inequality.  In its most elementary form, this inequality states that if $A = (a_{ij})$ is an $m\times n$ matrix with real entries, and we define the norm
$$\|A\| = \sup\biggl(\left|\sum_{i=1}^m \sum_{j=1}^n a_{ij}s_it_j\right| : |s_i| \le 1, |t_j| \le 1\biggr)$$
then
$$\biggl(\sum_{i,j} |a_{ij}|^{4/3}\biggr)^{3/4} \le \sqrt{2} \|A\|.$$
Are there more convincing examples of the importance of "exotic" values of $p$?  I remember wondering about this as an undergraduate but never pursued it.  As I think about it now, it does seem a bit odd from a pedagogical point of view that none of the textbooks I've seen give any applications involving specific values of $p$.  I didn't run into Littlewood's 4/3 inequality until later in life.
[Edit: Thanks for the many responses, which exceeded my expectations!  Perhaps I should have anticipated that this question would generate a big list; at any rate, I have added the big-list tag.  My choice of which answer to accept was necessarily somewhat arbitrary; all the top responses are excellent.]
 A: Hypercontractivity is a powerful technique that makes heavy use of $L^p$ spaces for $p \in (1,\infty)$.
Let $\|\cdot\|_p$ denote the $L^p$ norm.
Such results establish for an operator $T$ and function $f$ that
$$
\|f\|_q \leq \|Tf\|_p
$$
where $1<p < q$.
Perhaps the most striking example is the Bonami-Beckner inequality (originally due to Gross), which establish hypercontractivity for the Ornstein-Uhlenbeck operator and the noise operator, both parameterized by some variable $\epsilon$, on Boolean functions for appropriate values of $p, q$ and $\epsilon$.
The most famous application of the Bonami-Beckner inequality to the analysis of Boolean functions is the KKL inequality, which has had an enormous influence on the field.
Moreover, any time you see the words log-Sobolev inequality (which happens a lot when studying concentration of measure), hypercontractivity is lurking.
There are also reverse hypercontractivity results for $q < p < 1$. In particular, there are reverse versions for the noise operator and Ornstein-Uhlenbeck operator. These are used in the proof that the majority function is the most stable Boolean function (see here).
You can read more about hypercontractivity for Boolean functions in Ryan O'Donnell's book.
A: In fact $L^1$ and $L^\infty$ though natural in a naive sense are less well  behaved than the $L^p$ spaces for $1<p<\infty$ from the perspective
of dealing with PDEs.  Elliptic operators are not well behaved on the Sobolev spaces based on $L^1$ and $L^\infty$ while they are on the other $L^p$ spaces.
From the PDE perspective the more subtle friends of $L^1$ and $L^\infty$, namely
the Hardy space and its dual BMO of functions of bounded mean oscillation a have a good elliptic theory.  (See Stein's big book.)
A: Huge chunks of the theory of nonlinear PDEs rely critically on analysis in $L^p$-spaces.


*

*Let's take the 3D Navier-Stokes equations for example. Leray proved in 1933 existence of  a weak solution to the corresponding Cauchy problem with initial data from the space  $L^2(\mathbb  R^3)$. Unfortunately, it is still a major open problem whether the Leray weak solution is unique. But if one chooses the initial data from $L^3(\mathbb R^3)$, then Kato showed that there is  a unique strong solution to the Navier-Stokes equations (which is known to exist locally in time). $L^3$ is the "weakest" $L^p$-space of initial data which is known to give rise to unique solutions of the 3D Navier-Stokes.

*In some cases the structure of the equations suggests the choice of $L^p$ as the most natural space to work in. For instance, many equations stemming from non-Newtonian fluid dynamics and image processing involve the $p$-Laplacian 
$\nabla\left(|\nabla u|^{p-2}\nabla u\right)$ with $1 < p < \infty.$ Here the $L^p$-space
and $L^p$-based Sobolev spaces provide a natural framework to study well-posedness and regularity issues.



*

*Yet another example from harmonic analysis (which goes back to Paley and Zigmund, I think). Let
$$F(x,\omega)=\sum\limits_{n\in\mathbb Z^d} g_n(\omega)c_ne^{inx},\quad x\in \mathbb T^d,$$
where $g_n$ is a sequence of independent normalized Gaussians and $(c_n)$ is a  non-random element of $l^2(\mathbb Z^d)$. Then the function $F$ belongs almost surely to any $L^p(\mathbb T^d)$, $2\leq p <\infty$  and it does not belong almost surely to $L^{\infty}(\mathbb T^d)$.


There have been very recent applications of this resut to the existence of solutions to the nonlinear Schrodinger equations with random initial data (due to Burq, Gérard, Tzvetkov et al).
A: I feel as though this question may have come up before. Anyhow, the $\ell_4$ norm, and more generally the $\ell_{2k}$ norm for any positive integer $k$, come up naturally in Fourier analysis, since the $\ell_{2k}$ norm of the Fourier transform of $f$ equals the sum of $f(x_1)...f(x_k)\overline{f(y_1)...f(y_k)}$ over all $x_1+...+x_k=y_1+...+y_k$. That sort of sum comes up a lot in additive combinatorics, especially when $f$ is closely related to the characteristic function of a set. And you can get other norms by duality -- for instance the $4/3$ norm is the dual of the 4-norm, and therefore comes up too.
A: Tim, I've got two words for you: interpolation theorems (e.g., Riesz-Thorin and  Marcinkiewicz interpolation theorems). Such theorems let you pass from information about some operators on $L^1$ and $L^\infty$ to some operators on $L^2$ using all the intermediate exponents $p$.  
The point here is not that one actually cares about $L^{37.24}$ for its own sake, but the interpolation theorems show you that such "exotic" $L^p$-spaces can be at the service of her majesty $L^2$. I think for a student, these interpolation theorems provide an attractive reason to care about $L^p$ for all $p \geq 1$. 
This is not my area at all, so I welcome follow-up comments from analysts on this answer.
A: In PDEs, various values of p arise as degrees of regularity.  The Sobolev embedding theorems let you "trade in" generalized derivatives for classical derivatives.  You might need the exponent p to be above a certain threshold to get a desired regularity result.
Still, I agree with your observation that much of the time the values of p that matter are 1, 2, and infinity.
A: In probability, the $L^p$ norms give you the $p$-th moments of a random variable, and the relationship between them can tell you a lot about its distribution.  For example, the 3rd and 4th moments tell you something about how symmetric and how concentrated about its mean a distribution is.  Statisticians give them cool names like "skewness" and "kurtosis".
I'll also mention a recent startling theorem of Nualart et al, which says that a sequence of random variables taken from a Wiener chaos converge in distribution to a certain limit if and only if their fourth moments are converging to the right thing.  (First and second moments are not sufficient.)
A: $L^p$ norms for $p$ large (but $<\infty$) have played a crucial role in additive combinatorics via the technique of almost-periodicity (most recently in the spectacular bounds for three-term arithmetic progressions by Kelley and Meka).
The reason is the following: for many problems in additive combinatorics we want to understand how some object like $\langle 1_A, f\rangle$ changes if we translate $A$ by some $t\in G$, where $G$ is a finite abelian group and $A\subseteq G$ has density $\alpha=\lvert A\rvert/\lvert G\rvert$.
This means we want to bound
$$ \lvert \langle 1_{A-t},f\rangle - \langle 1_A,f\rangle\rvert = \lvert \langle 1_A, \tau_t f-f\rangle\rvert,$$
where $\tau_tf(x)=f(x+t)$. By Hölder's inequality this is at most, for any $1\leq p\leq \infty$,
$$ \| 1_A\|_{p/(p-1)}\| \tau_t f-f\|_p=\alpha^{1-1/p}\| \tau_t f-f\|_p.$$
The natural scale for the left-hand side is $\approx \| 1_A\|_1\|f\|_\infty=\alpha \| f\|_\infty$, and say we only want to save some $\epsilon$ over this trivial bound, so our goal is to bound the right-hand side by $\ll \epsilon \alpha \|f\|_\infty$, say.
Almost-periodicity techniques (originating in work of Croot and Sisask) allow us to bound (not generally, but for many $f$ of interest) $\| \tau_t f-f\|_p\ll \epsilon \|f\|_\infty$ for 'many' $t$, where 'many' is some explicit quantity with good (polynomial) dependence on $p$ and $\epsilon$.
Now if one just tries this with $p=2$ then you need to save some factor of $\alpha^{-1/2}$ somehow. But the magical trick is if we choose $p\approx \log(1/\alpha)$ then $\alpha^{-1/p}\ll 1$, and hence one can find 'many' $t$ such that $\langle 1_{A-t},f\rangle=\langle 1_A,f\rangle+O(\epsilon \alpha \|f\|_\infty)$ where 'many' depends polynomially on $\epsilon$ but only logarithmically on $\alpha$ (previous approaches via Fourier analysis just used the $p=2$ norm and had polynomial dependence on $\alpha$ as a result).
(This is very similar to the use of large $p$ in Banach space theory as mentioned in Bill Johnson's answer.)
A: Tim, here is one very specific example that a computer scientist who cares only about $L_1$ and $L_2$ should find appealing. The norm of  $\ell_1^n$ is, up to a constant, the same as the $\ell_p^n$ norm when the conjugate index to $p$ is $\log n$.  As was already mentioned in this thread, the $L_p^n$ norm is uniformly convex when $1<p<\infty$ and the modulus of convexity is known.  This fact, frequently used by researchers in Banach space theory, was used by Lee and Naor to give a strikingly simple proof of the Brinkman-Charikar result on the impossibility of dimension reduction in $L_1$.
A: Although nonlinear PDE's are mentioned by John Cook, he seems to still concede that $p = 1, 2, \infty$ are the most important. I beg to differ. I will give only one specific example that I am familiar with. As others have noted, Terry Tao has written a lot about using $L_p$ estimates to study other types of nonlinear PDE's.
In the 70's and 80's, there were major breakthroughs in the use of elliptic PDE's to prove global theorems in geometry and topology by Sacks-Uhlenbeck and Schoen-Simon-Yau in minimal surfaces, Uhlenbeck, Taubes, and Donaldson in Yang-Mills theory, and Gao and Anderson-Cheeger in Einstein manifolds. The critical tool used here were sharp Sobolev inequalities and $L_p$ estimates of the first derivative to the solution of a nonlinear PDE. The technical tool that is often used is called Moser iteration, where initial $L_p$ bounds on the gradient are bootstrapped into stronger bounds on the solution. These types of estimates can also be applied to nonlinear parabolic PDE's, including the Ricci flow.
All of this has led to a tremendous growth in the study of nonlinear elliptic PDE's, along with their applications to global differential geometry and topology, as well as mathematical physics. The $L_p$ theory, where $p \ne 2$, plays a crucial role in most of this work.
A: I think a nice example are Lieb--Thirring inequalities. Consider the Schroedinger operator $H = \Delta + V$ with potential $V \in L^{\gamma + d/2}(\mathbb{R}^{d})$, where $d \geq 3$. Then H defines an (unbounded) operator on $L^2(\mathbb{R}^d)$, whose essential spectrum is $[0,\infty)$ and which has negative eigenvalues $E_j$ (countably many). The Lieb--Thirring inequality then tells us
$$
 \sum |E_j|^{\gamma} \leq const  \|V\|_{L^{\gamma + d/2}(\mathbb{R}^{d})}.
$$
This inequalities requires $L^p$ for $p \in (0,\infty)$.
There are other examples, but they are somewhat more technical to state ...
A: One answer that I don't see explicitly mentioned here is that $L^\infty = \lim_{p\to\infty} L^p$, so to speak.  So if you in fact care about $L^\infty$ but there are easier theorems regarding $p<\infty$, then do what we do in analysis, and get a sequential approximation.
A: Here is an algebraic answer to your question.
See a related answer of mine for more details.
First, as I explain in the answer cited above, it is very natural
to replace the number p by its reciprocal 1/p, i.e., define L_p := L^{1/p}.
Thus L^1, L^2, L^∞ are denoted by L_1, L_{1/2}, and L_0 in this notation.
For an arbitrary measurable space Z (i.e., a commutative von Neumann algebra),
and, more generally, for an arbitrary noncommutative measurable space Z (i.e.,
a noncommutative von Neumann algebra) we can define the space L_p(Z) for all p∈CP,
where CP is the set of all complex numbers with a nonnegative real part.
Note that no choice of measure (weight) on Z is needed to construct L_p(Z).
Note that L_0 just consists of bounded functions on Z and L_1 consists of finite
complex-valued measures (weights) on Z.
These spaces are graded components of a complex unital CP-graded *-algebra.
In a certain precise sense one can say that this CP-graded *-algebra is a free algebra
generated by all bounded functions on Z in grading 0 and all finite complex-valued measures (weights) on Z in grading 1,
with the obvious relations coming from the Radon-Nikodym theorem
and (in the noncommutative case) the modular automorphism group
(which essentially explains how weights (i.e., noncommutative measures)
commute with bounded functions).
If ℜp∈[0,1] then L_p(Z) is a Banach space, otherwise it is a quasi-Banach space.
Also, if ℜp∈[0,1], then L_p(Z) can be obtained as the complex interpolation of L_0(Z) and L_1(Z)
corresponding to the parameter p.
A lot of theorems in the (noncommutative) integration theory can be proved by simple
algebraic manipulations in this CP-graded *-algebra.
These algebraic manipulations require that one has access to all graded components,
not just components with gradings 0, 1/2, and 1.
Let me give just one example for L_p spaces where p is imaginary.
Suppose μ is a weight on M
(if μ is bounded, i.e., μ(1)<∞, then μ∈L_1(Z), otherwise we should think of μ as an element
of the extended positive cone EL_1^+(Z)).
Suppose furthermore that t is an imaginary number and x∈L_0(Z).
Then we have μ^t∈L_t(Z), x∈L_0(Z), μ^{-t}∈L_{-t}(Z)
and their product is σ^μ_t(x) := μ^t x μ^{-t} ∈ L_0(Z).
The one-parameter automorphism group
σ^μ is called the modular automorphism group of the weight μ
and it is a very important notion of noncommutative geometry.
(In the commutative case we always have σ^μ_t(x)=x.)
A: The last paragraph of section 6.1, "Basic theory of $L^p$ spaces", in Folland's Real Analysis neatly summarizes a lot of the points made in other answers:

We conclude this section with a few remarks about the significance of the $L^p$ spaces. The three most obviously important ones are $L^1$, $L^2$, and $L^\infty$. With $L^1$ we are already familiar [from the development of Lebesgue integration in earlier chapters]; $L^2$ is special because it is a Hilbert space; and the topology on $L^\infty$ is closely related to the topology of uniform convergence. Unfortunately, $L^1$ and $L^\infty$ are pathological in many respects, and it is more fruitful to deal with the intermediate $L^p$ spaces. One manifestation of this is the duality theory in section 6.2; another is the fact that many operators of interest in Fourier analysis and differential equations are bounded on $L^p$ for $1 < p < \infty$ but not on $L^1$ or $L^\infty$.

A: I had a similar question when I was first learning about the Lebesgue spaces: does anyone actually use these spaces when p<1?  There are obvious technical problems with these spaces since the unit balls are not convex; despite this fact, the answer is yes.  There are a number of interesting multilinear operators which are bounded maps from $L^{p_1}\times L^{p_2}\times ...\times L^{p_n}\rightarrow L^r$ where the exponents satisfy the condition $\displaystyle{\sum_{j=1}^n\frac{1}{p_j}=\frac{1}{r}}$.  Now, if $n\ge 3$ and $p_i=2$ for all i, then this forces r to be some fraction less than 1.  So, even if one only cared about $L^2$ and $L^\infty$, one would still run into values of $p<1$.
One such class of operators is comprised of multilinear variants of the Hardy-Littlewood maximal operator.  They were studied fairly recently by Ciprian Demeter, Terence Tao, and Christoph Thiele, e.g. in Maximal multilinear operators, Trans. Amer. Math. Soc. 360 (2008), 4989-5042, doi:10.1090/S0002-9947-08-04474-7 arXiv:math/0510581.  Another type of operator in this spirit is the Biest operator studied by Camil Muscalu, Terence Tao, and Christoph Thiele (e.g. $L^p$ estimates for the biest II. The Fourier case, Math. Ann. 329 (2004) 427–461, doi:10.1007/s00208-003-0508-8, arXiv:math/0102084).
A: I guess that $L^1$, $L^2$, and $L^\infty$ just seem natural because they are so intimately related to obvious everyday concepts -- sums, averages, maxima, root-mean-square (Hilbert space, ...).  So I doubt that the occasional theorem that involves another $L^p$ space explicitly will ever convince anyone that any other $L^p$ space is equally significant.
But let me just mention another one of them anyway, a marvelous theorem of Beurling: Let $p\in(1,\infty)$.  Then the family of functions of the form $f(x) = \sum_{k=1}^n a_k\rho(\theta_k/x)$---where $\rho(x) = x -\lfloor x\rfloor$ and $\sum a_k\theta_k = 0$---is dense in $L^p(0,1)$ iff the Riemann Zeta function has no zeros in the half-plane $\sigma > 1/p$.
A: $L^0$, which for a finite measure is just the set of all measurable functions, is important in probability.  When you're modeling some physical phenomenom, there is often no canonical choice of probability measure, so you sometimes work with a class of equivalent probabilities, i.e. ones that have the same null sets.  The only two $L^p$ spaces that are invariant under changing to an equivalent probability are $L^\infty$ and $L^0$.  The former space is often too small for modeling purposes, and so you are forced to work with $L^0$.  It's topologized by convergence in measure/probability, and it's pretty horribly non-convex.
A: Here is an example from probability theory. Let $X_i$ be a sequence of independent identically distributed random variables in $L^1$. The strong law of large numbers asserts that the mean converges to the expectation.
$$a.e. \quad {1\over n}\ \sum_{k=0}^{n-1} X_k \rightarrow E(X_0).$$
What can be said about the speed of convergence ? If we assume that the $X_i$ are in $L^p$ for some $p\in ]1,2[$, then we have:
$$a.e. \quad {1\over n}\ \sum_{k=0}^{n-1} X_k = E(X_0) + o(n^{1/p-1}).$$
