Introduction
The OP mentioned science majors, and since there are loads of examples from science, we might as well start (and probably finish) there. There is a mathematical viewpoint about the utility of bases and abstract vector spaces, but I'll leave that aside and stick to practical concerns of scientific computing. Real examples are complicated and not just by their detail. Pedagogically they are complicated by drawing on bits of linear algebra from throughout a standard first course and sometimes from extra bits you pick up as needed in your work after the standard course. They are further complicated by being based on principles of science, each of which might be unfamiliar to some of one's students. I think one can get around these complicating factors and motivate students. If the problem, including the science and the basic role of linear algebra, is easy enough to understand, then the problem might tempt their curiosity. The main point of these examples will be that when basis vectors are meaningful, their need and usefulness seems unquestionable. In the last group of examples, they are less meaningful but their utility is apparent; and it is less easy to see the vectors, which are functions, as $n$-tuples in a natural way.
$\def\y{{\bf y}}\def\x{{\bf x}}\def\A{{\Bbb A}}\def\R{{\Bbb R}}$
The problem with $\R^n$ that the OP describes has several sources.
- The vectors and their coordinates are indistinguishable.
- Students have been dealing with $\R^2$ for years.
- The standard basis is orthonormal, which is a most convenient sort of basis.
Any other basis looks pretty inconvenient to a beginner. So $\R^n$ feels "canonical." Furthermore:
- A main purpose of a basis is to make computation possible by mapping a vector space to $\R^n$.
- To a practically minded person, the real work, the computation, is done in $\R^n$, and that is where the real importance lies.
- So why not do this in the easiest, most obvious way possible? For instance, use the coefficients of polynomials.
- It would help the student to accept a basis if the basis vectors had a meaning or represented elementary concepts.
- The standard basis in $\R^n$ for $n=2,3$ has such a meaning: They tell you how to plot a point in a most simple way. And that meaning is ingrained in the beginner's mindset.
One extra thing I would point out to the students: You won't be doing much number-crunching yourself. You'll be setting up computations to solve problems. It will be easier to set the problem up in a convenient basis and let the computer figure out how to get back to $\R^n$.
In the scientific examples below, I've selected examples in which a basis vector corresponds to a fundamental element or structure in the problem being modeled. In function approximation, in which way too "many bases appear," this is less the case. But these are used in important applications. For instance, in radial basis function interpolation (already a mouthful), "surface reconstruction" may sound abstruse, and students might need help imagining what is so important about modeling a surface. Well, monitoring the growth of a cancer tumor or visualizing the surfaces in a heart seem important applications to me.
Linear Approximations
Science, to speak broadly, uses mathematics to model phenomena, often under simplifying assumptions that make the model approximate. The mathematics itself often needs to be computed approximately. Linear approximations are the simplest and are remarkably robust and accurate for solving real problems.
Spectra (chemistry, physics, group theory)
In a comment I made soon after the question was asked (now deleted), I mentioned that the normal/fundamental modes form a basis for the small vibrations of a molecule or other oscillatory system (small = linearized approximation). In the linearized approximation, the oscillations are given by (complicated) functions of time representing the displacements of coordinates from the rest position. These oscillations form a linear system, a vector space. The remarkable thing is that any vibrational motion is a linear combination of the basis of normal modes, each of which may have its own frequency. This applies to spring-mass systems (@Mark H) as well.
Normal-mode decompositions are used in spectroscopy in chemistry, say, to identify compounds in a mixture. This use of spectroscopy will have been seen by all the chemistry majors at my college by the end of their sophomore year, although they will have had little explanation of the mathematics. They'll see it again in upper-level chemistry along with group theory, where they might learn this gem that was a highlight of my undergraduate education: The irreducible representations of the vibrational symmetries of methane (or other tetrahedral molecule) cannot be computed without the complex numbers. That is to say, the real phenomenon cannot be fully understood with just real numbers. As a liberal arts major, I could rejoice that what was real was complex.
Something I learned later was the significance of the two complex-conjugate irreducible representations that appear in the analysis of methane. The vibrational modes associated to these symmetries are the only ones that generate an oscillating dipole that is necessary for infrared absorption. There are three modes associated to each of the two symmetries, both labeled $t_2$ below but with different frequencies. Further, they oscillate at infrared frequencies (wavenumbers $3156.8\ \text{cm}^{-1},\ 1367.4\ \text{cm}^{-1}$. It is why methane is a powerful greenhouse gas.
Approximating matrices
When you have data, such as the oscillations of the coordinates of the atoms in a molecule as in the preceding example -- the $\R^n$ representation -- the standard coordinates don't necessarily in themselves represent an important concept in the problem, unlike the normal-mode coordinates, which indicate how much of each fundamental element, a normal mode, makes up a given vibration. In the vibration problem, there are as many normal modes as degrees of freedom, and each normal-mode vector is equally important. In some problems, some vectors are more important than others, and identifying them can reveal the structure of phenomena. We'll look at one way to do that.
It turns out that any matrix $A$ can be written
$$A = U S V^*\,,$$
where
- $V^*$ is the transpose of $V$ (or conjugate transpose if $A$ is complex),
- $A$ is an $m\times n$ matrix,
- $U$ and $V$ are $m \times m$ and $n \times n$ orthogonal matrices (or, resp., unitary ones if $A$ is complex -- that is, in both cases, $V^*=V^{-1}$),
- and $S$ is a $m \times n$ diagonal matrix with nonnegative real diagonal entries
$$s_{1} \ge s_{2} \ge \cdots \ge s_{n} \ge 0 \,.$$
If $r$ is the rank of the matrix $A$ and $r < n$, then we have
$$s_{1} \ge s_{2} \ge \cdots \ge s_{r} > 0 \ \text{and}\ s_{r+1}=s_{r+2}=\cdots=s_{n}=0\,.$$
If the real numbers
$$s_{1}, s_{2}, s_{3},\dots$$
decrease rapidly, we can replace the small ones by zero and change the matrix $A$ by only a small amount.
This factorization of $A$ is called the singular value decomposition or SVD and the nonzero $s_j$ are called the "singular values" of $A$. The SVD factorization is not always covered in a first linear algebra course. It could be stated as a fact and used with the examples below to motivate future learning. It is worth considering simply because it's a really cool structure theorem that is in widespread use to solve real problems. Have the students google "principal component analysis" or "PCA."
The SVD is readily available in any decent system or library for linear algebra, and students may have to use it without ever covering it in class:
Python,
R (svd()
),
Sage,
Maple,
Matlab,
Mathematica,
Maxima
GNU GSL,
etc.
We'll give a brief explanation of how it identifies "principal" components of a linear model.
First, let's break down the factorization of $A$.
The matrix $A$ represents a linear mapping from $\R^n$ to $\R^m$.
The columns of $U$ and $V$ are orthonormal bases (the best kind) for the output $\R^m$ and the input $\R^n$, respectively. First some notation: for a matrix $M$, let $M_{ij}$ denote the $ij$ entry (two indices) and $M_j$ denote $j$-th column matrix (one index). Then to compute $Ax$ for a vector $x$, we may write
$$Ax = \Big(\sum_{j,k} U_{ij} s_j V_{kj}x_k\Big)_{i=1}^m = \sum_j (s_j\; (V_j^* \cdot x)) U_{j} \,. \tag{1}$$
Thus we first compute $V^*x$, which are the coordinates of $x$ with respect to the basis given by the columns of $V$. Then $SV^*x$ scales these coordinates by the $s_j$. These scaled coordinates are the coordinates of $Ax$ with respect to the basis given by the columns of $U$, that is $Ax = U(SV^*x)$. It's pretty amazing that every matrix action can be viewed in terms of choosing two orthonormal bases and scaling coordinates. We can also see what happens when we replace $s_j$ by $0$ for a single value of $j$. It is in effect the same as deleting the $j$th term from (1), which introduces an error bounded by
$$\|(s_j\; (V_j^* \cdot x)) U_{j}\| \le s_j\,\| x\| \,,$$
since $\|U_j\|=\|V_j\|=1$ and $s_j \ge 0$.
If we replace all $s_j$ by $0$ for $j >J$, then the error relative to $\|x\|$ will be bounded by
$$\sum_{j >J} s_j \,. \tag{2}$$
Image processing (data compression)
I'm counting it as "science" because we'll use the image of a scientist. Mainly, I introduce the notion this way because it's visual, and it's hard to believe one's eyes. In this application, we're approximating a matrix itself. Since the matrix entries are extracted from image data and vice versa, students can see the approximation improve. The initial matrix entries are the grayscale values (between $0$ and $1$) of the grid of pixels. The singular values this matrix decrease from $s_1=338$ to $s_4=30$, to $s_{50}=3.3$, to $s_{100}=1.4$, and finally to $s_{841}=0.039$ -- a decrease of five orders of magnitude.
Below are the images from the matrix
$$A^*=\sum_{j=1} ^J s_j\, U_{j} (V_j)^T$$
for the singular values up to $J=4,\,50,\,100,\,841$ (all). The relative error (RE) was calculated using matrix 2-norms, $\text{RE} = \|A^*-A\|/\|A\|$, but one should expect a weak correlation between the error of approximating a linear map and how accurate the picture looks. The Frobenius norm gives relative errors of 21%, 5.6%, 3.6%, and 0% respectively, which perhaps seem closer to how the images appear. All norms give estimates much lower than the upper bound $\sum_{j>J} s_j$ deduced above.

Here are a couple of codes for constructing the image approximation (I used this image from Wikipedia, scaled by a factor of 1/3 and converted to grayscale):
(* Mathematica *)
img = Import["<path/url to image file>"]; (* convert to grayscale if necessary *)
{U, S, V} = SingularValueDecomposition[ImageData@img]; (* assuming "GrayScale" image *)
J = 4; (* change as desired *)
Image[U[[All, ;; J]] . S[[;; J, ;; J]] . Transpose@V[[All, ;; J]]]
% MATLAB
img = imread('<path to image file>'); % convert grayscale if necessary
[U,S,V] = svd(im2double(img)); % convert grayscale img to double
J = 50; % change as desired
newimg = U(:,1:J) * S(1:J,1:J) * V(:,1:J)';
imshow(newimg)
DNA analysis (biology)
Quantitative Understanding in Biology: Principal Component Analysis (lecture notes by Ju et al.) is a wonderful introduction to principal component analysis (PCA), which may be read by any student in linear algebra. It refers the connection to the singular value decomposition to a technical paper, but it does a good job explaining why identifying principal components and using them as basis vectors helps one tame an unwieldy dataset. The authors describe how some common DNA datasets have 10,000+ measurements per sample and others around 600,000-900,000. They point out, "These datasets are typically never fully visualized because they contain many more datapoints than you have pixels on your monitor." They introduce the notion of PCA with a simple numerical example, consisting of a dataset of two samples resulting from the measurement of 5 genes ($g_1$ through $g_5$), in which there is hidden correlation. They show that the 5-dimensional representation may be reduced to 2 dimensions, using the principal components. Further they explain the weakness of the standard coordinates: "knowing the expression level of $g_1$ will not give us any information about other genes." And they explain the significance of the principal components to understanding the biology:
However, in reality the expression levels of multiple genes tend to be correlated to each other (for example, pathway activation that bumps up the expression levels of $g_1$ and $g_2$ together, or a feedback interaction where a high level of $g_3$ suppresses the expression level of $g_4$ and $g_5$), and we don't have to focus on the expression levels individually. This would mean that by knowing the expression level of $g_1$ we can get some sense of the expression level of $g_2$, and from the level of $g_3$ we can guess the levels of $g_4$ and $g_5$.
Thus the principal basis vectors reflect biological function. That's pretty cool.
Protein Conformational Motions of HIV-1 Protease (biology, medicine, biomedical)
In Singular Value Decomposition of Protein Conformational Motions: Application to HIV-1 Protease, Teodoro et al. analyze the conformational changes of an enzyme via the SVD. "An illustrative example of these changes is the opening and closing of the binding site of HIV-1 protease (HIV Pr)...." HIV Pr plays an important role in the life cycle of the HIV virus, and the opening and closing of the binding site allows the enzyme to do its job. Disrupting its function is the focus of protease inhibitors. The analysis of the motions is similar to the molecular spectra discussed above, except proteins exhibit both small vibrations and larger motions, and the interest lies in the larger motions. The dimension of the space of motions is $597$, and the SVD was used to focus on a few principal basis vectors.
The left singular vectors [i.e., the columns of $U$] obtained from the decomposition correspond to modes of collective motion whose displacements are directly proportional to the value of the corresponding singular values. The calculated right singular vectors [i.e., the columns of $V$] correspond to the projection of the original trajectory in the new basis.
In this way, they were able to show that "the protein jumps rapidly between stable intermediate conformations until it reaches the open conformation." Again, the problem cannot easily be investigated in terms of $597$-tuples, and choosing the right basis, in which the vectors reflect structures within the problem, allows one to deduce how those structures work.
Function approximation
Function approximation is not a science topic per se, but it is used broadly in numerical solutions to scientific problems. Let's say you have a bunch of function values $y^{(j)}$ at a number of points $\x^{(j)}=(x_1^{(j)},x_2^{(j)},\dots,x_n^{(j)})$ for $j=1,2,\dots,m$, and you'd like to construct a continuous function $f(\x)=f(x_1,x_2,\dots)$ based on this data. How to do it?
Viewed as a linear problem, we would choose a basis of functions $\{\phi_1,\dots,\phi_n\}$
and solve for coefficients $c_j$ such that
$$f(\x) = \sum c_j \phi_j(\x)$$
meet some criteria. The interpolation problem is to find a function that passes through each point so that $f(\x^{(j)})=y^{(j)}$. The fitting problem is to find a function that minimizes the error between the values $f(\x^{(j)})$ of the function and the values $y^{(j)}$ at each point; the error typically used is least squares, which problem can be solved with linear algebra. When the values $y^{(j)}$ are the values of a given function, then the problem is one of approximating the given function. If the basis functions have the property that $\phi_k(\x^{(k)})=1$ and $\phi_k(\x^{(j)})=0$ for $j\ne k$, then the basis is particularly nice for the interpolation problem: $c_j=y^{(j)}$ and there's no need to solve a system of equations for $c_j$. Such a basis is called a cardinal basis.
Polynomial interpolation
The power basis $1,\ x,\ x^2,\dots, x^n$ is often numerically poorly conditioned for the interpolation problem. The Chebyshev polynomials $T_j$ form a better basis, and the Lagrange basis, which is a cardinal basis, is better still in its barycentric form (see Berrut & Trefethen, eq. (4.2)). There are other families of orthogonal polynomials besides Chebyshev that have their uses.
Chebyshev approximation
A truncated power series (Taylor polynomial) is often not a good approximation away from its center of expansion. However, for a function that is analytic over a finite interval $[a,b]$, its Chebyshev series converges rapidly on $[a,b]$ with an error bound that is easy to state and apply practically. Let $f(x)$ have a Chebyshev series
$$f(x) = \sum_{j=0}^\infty c_j T_j\left( \frac{2 x - (a+b)}{b-ab} \right) \,.$$
If the series is truncated at the term $j=J$, then since $|T_j| \le 1$ on the interval, the truncation error is bounded by
$$\sum_{j>J} |c_j| \,, \tag{3}$$
which is similar to the SVD truncation error bound in (2). If the coefficients $c_j$ are rapidly decreasing,
then $c_J$ is usually a good bound on the error.
If we interpolate $f(x)$ at the Chebyshev points $x_k = \cos(\pi k /J)$, $k=0,1,\dots,J$
using the Chebyshev basis, we approximate the truncated Chebyshev series
$$p(x) = \sum_{j=0}^J a_j T_j\left( \frac{2 x - (a+b)}{b-ab} \right)
\approx \sum_{j=0}^J c_j T_j\left( \frac{2 x - (a+b)}{b-ab} \right)
\approx f(x)\,.$$
The Chebyshev basis is an orthogonal basis with respect to the inner product $\langle f, g \rangle = \int_a^b f(x) g(x) / \sqrt{(x-a)(b-x)}\, dx$, and truncating the series is equivalent to orthogonal projection onto
a finite dimensional subspace (= polynomials of degree at most $J$). As such, the truncated series is an optimal solution to approximating the function $f(x)$ with a polynomial of degree at most $J$ with respect to the norm induced by the inner product.
Further it is easy to evaluate, integrate, differentiate and find the roots of Chebyshev interpolations.
It is one of the best ways to solve transcendental equations (see Boyd (2014)).
Chebyshev approximation is the basis of the Chebfun project.
Fourier approximation
Virtually the same thing can be said for Fourier approximation, which in fact might be viewed as the basis for Chebyshev approximation. The Fourier basis, mentioned by @Kuba, for the interval $[-\pi,\pi]$ is $\phi_j(x) = e^{i j x}$, $j\in{\Bbb Z}$. Since $|e^{i j x}|=1$ and the coefficients of Fourier series for analytic functions are rapidly decreasing, we get similar approximation properties, including orthonormality of the basis (with respect to the inner product $(2\pi)^{-1}\int_{-\pi}^\pi f(x) g(x) \, dx$). One reason this and the Chebyshev approximation are important is that the Fast Fourier Transform (FFT) makes these approximations computationally efficient.
Sinc interpolation (engineering, signal processing, image processing)
The sinc function $(\sin t)/t$ or the normalized sinc function $\sin(\pi t)/(\pi t)$ (for $t\ne0$) can be used to form a cardinal basis for interpolating sampled signals or other regularly spaced time series. It should be understood that the sinc functions are taken to be continuous and have the value $1$ at $t=0$.
Suppose we have a regular sample rate of $1/h$ and thus the sampled times are $t_j=jh$.
The sinc basis functions are
$$\phi_j(jh) = 1\,,\ \text{otherwise}\ \phi_j(t) = {\sin(\pi (t-jh)) \over \pi (t-jh)} \,.$$
The interpolation formula for the sampled data $y^{(j)}=f(jh)$ of a signal $f(t)$ is
$$f(t) \approx \sum_j y^{(j)} \phi_j(t) \,.$$
Radial basis function interpolation (neural networks, engineering, signal processing)
Radial basis function (RBF) interpolation can be used to interpolate "unstructured" data of any dimension (the sampling nodes need not be on a regular grid). For instance, the sampling nodes could be the locations in a network of weather stations or other geographic phenomena. Applications include the approximation of functions on arbitrary surfaces; surface reconstruction from medical imaging or geological stratigraphic data; machine learning; and so forth.
An RBF approximation is built on a radial basis $\phi_j(x)=\phi(\| \x-\x_j\|)$ and has the form
$$f(\x) \approx \sum_{j=1}^J c_j \phi(\| \x-\x_j\|) \tag{4}$$
for some function $\phi(r)$ of the distance from a sampling node, typically a Gaussian $e^{-\alpha^2 r^2}$ for some "shape parameter" $\alpha$.
Here $\x$ may be in $\R$ or $\R^n$.
Solving for the coefficients $c_j$ is a linear problem and might be done either as a least-squares fitting problem or as an interpolation problem.
As an example,
Maz'ya & Schmidt (2007)
and
Boyd & Wang (2009)
independently derived a cardinal basis for Gaussian RBF approximation on an infinite grid with uniform spacing $h = x_j - x_{j-1}$ ($-\infty<j<\infty$)
and shape parameter $\alpha$:
$$\phi_j(x)=\phi(x-x_j),\quad\text{where}\quad
\phi(x)= {\alpha ^2\over\pi}\,
{\sin \left({\pi x}/{h}\right) \over {\sinh}\left({\alpha^2 x}/{h}\right)} \,.$$
Since for this grid, the sum would be infinite, one usually truncates it and sums over a finite interval.
(It is related to sinc interpolation but the $x$s and $h$s in $\sin \left({\pi x}/{h}\right)/(\pi x/h)$ and ${\sinh}\left({\alpha^2 x}/{h}\right)/(\alpha^2 x/h)$ cancel out.)