Why does resonance occur at only standing wave frequencies in a fixed string?

One way to explain it is to look at the string as a waveguide. You can send wave packets over the string, and they're reflected at the ends:

Now, you can send in wave packets repeatedly rather than just once:

This is essentially exciting the string with a (non-sinuoidal) periodic signal of pulses. But it is not resonance: the new pulses don't amplify the old ones in a meaningful way, they just randomly add or cancel.

However, when you send the next pulse always after exactly one period of the string, it will consistently add to the old pulse's amplitude, and that's where you get resonance: a small input amplifies to a large oscillation over time.

This works ^{_{(in case of a wave equation with linear dispersion)}} for any periodic signal, i.e. for a regular stream of pulses as well as a sine input:

...Again, only if you hit the correct frequency. If not, the behavious is something like that:

Source code: (Haskell with dynamic-plot)

{-# LANGUAGE TypeFamilies #-}

import Graphics.Dynamic.Plot.R2
import Text.Printf

import Data.VectorSpace
import Data.VectorSpace.Free

import qualified Data.Vector.Unboxed as VU

import Data.Semigroup

rk₄ :: (VectorSpace v, t ~ Scalar v, Fractional t)
            => t -> (t -> v -> v) -> v -> [(t,v)]
rk₄ h f y₀ = go 0 y₀
 where go t y = (t, y) : go (t+h)
          (y ^+^ (h/6)*^(k₁ ^+^ 2*^(k₂^+^k₃) ^+^ k₄))
        where k₁ = f  t        y
              k₂ = f (t+h/2) $ y^+^(h/2)*^k₁
              k₃ = f (t+h/2) $ y^+^(h/2)*^k₂
              k₄ = f (t+h)   $ y^+^h*^k₃

type Time = Double
type Amplitude = Double
type Velo = Double

newtype StringState
  = StringState { getFDSState :: VU.Vector (Amplitude, Velo) }

instance AdditiveGroup StringState where
  zeroV = StringState $ VU.replicate 128 zeroV
  negateV (StringState ys) = StringState $ VU.map negateV ys
  StringState ys ^+^ StringState zs
      = StringState $ VU.zipWith (^+^) ys zs

instance VectorSpace StringState where
  type Scalar StringState = Double
  μ *^ StringState ys = StringState $ VU.map (μ*^) ys

spatDeriv² :: StringState -> VU.Vector Double
spatDeriv² (StringState ys)
 = VU.singleton 0
   <> VU.zipWith3 (\(l,_) (m,_) (r,_)
                    -> (l + r - 2*m)*fromIntegral (VU.length ys)^2)
             ys (VU.drop 1 ys) (VU.drop 2 ys)
   <> VU.singleton 0

stringDynamics :: (Time -> Amplitude) -> StringState
                     -> [(Time,StringState)]
stringDynamics excite y₀
   = sparsen $ rk₄ (1/(20*fromIntegral oversampling)) propag y₀
 where propag t st@(StringState ys) = StringState
           $ VU.zip
              (VU.map snd ys
                 VU.// [(0, 9 * (excite t - fst (ys VU.! 0)))])
              (VU.map (^*(0.1)) (spatDeriv² st))
       oversampling = 4
       sparsen (q:qs) = q : sparsen (drop (oversampling - 1) qs)

stringPlot :: StringState -> DynamicPlottable
stringPlot (StringState ys)
  = lineSegPlot $ zip [0, 1/fromIntegral (VU.length ys)..]
                      (map fst $ VU.toList ys)

period :: Time
period = 2*sqrt 10

animateResponse :: (Time -> Amplitude) -> DynamicPlottable
animateResponse excite
  = plotLatest
     [ legendName (printf "t = %.1g" t) $ stringPlot ys
     | (t,ys) <- stringDynamics excite zeroV
     ]

main :: IO ()
main = do
 plotWindow [
   animateResponse $ \t -> let ccl = floor $ t/period
                               t' = t - fromIntegral ccl*period
                           in exp (-8*(t'-5)^2)/3
  ]
 plotWindow [
   animateResponse $ \t -> sin (2*pi*t/period)/4
  ]
 plotWindow [
   animateResponse $ \t -> sin (2*pi*t*1.20124)/4
  ]
 return ()

Are you familiar with the steady-state solutions to the equation $$ \frac{d^2 x}{dt^2}+\Omega_0^2 x= F \sin (\Omega t)? $$ (If not, look for a solution $x(t) \propto \sin (\Omega t)$ and find the constant of proportionality)

When you solve this you can immediately see that you can excite an oscillator with frequency $\Omega$ even when $\Omega$ is not equal to the resonant frequency $\Omega_0$. You can then apply this to each mode in the normal mode decomposition of the string and see that the string will vibrate at the applied frequency $\Omega$, but not very much unless $\Omega\approx \Omega_n$ where $\Omega_n$, for some $n$, is the frequency of the $n$-th mode.

Note added: To expand on what I said above: if you apply a force $F(x) \sin \omega t$ you first decompose into normal modes $$ y(x,t)= \sum_n a_n(t) \sin\left (\frac{n\pi x}{L}\right) $$ and then each mode $a_n(t)$ obeys $$ \frac{d^2 a_n(t)}{dt^2}+\Omega_n^2 a_n(t)= F_n \sin (\Omega t) $$ where $\Omega_n$ is the frequency of the mode $n$ and $ F_n$ is the $n$-th Fourier coefficient of $F(x)$.