arithmetic operations don't check compatibility of spectral_axis

Arithmetic operations succeed so long as the flux arrays are of compatible shape, but do not check for equivalent spectral axis values. For example,

```
from specutils import Spectrum1D
import numpy as np
from astropy import units as u

sp1 = Spectrum1D(np.random.random(10)*u.Jy, spectral_axis=np.linspace(0,1,10)*u.pix)
sp2 = Spectrum1D(np.random.random(10)*u.Jy, spectral_axis=np.linspace(2,8,10)*u.pix)
sp2 - sp1
```
succeeds but returns a `Spectrum1D` object with a spectral axis adopted directly from `sp2`. Instead this should raise an error (or perhaps eventually support optional interpolation for cases in which there is sufficient overlap).

Related: #198

Start by running the Spectrum1D example in the issue and tracing the arithmetic operation that currently accepts different spectral_axis values. Done means incompatible spectral axes are rejected consistently; the issue leaves optional interpolation and the relevant files or tests unspecified.

Written by the indexing model from the issue text.