Q:

Dissipative time-stepping scheme for first order in time system

When solving semi-discrete equations (originating from finite element models, for example), which are second-order in time of the form
\begin{equation}
M\ddot d + C\dot d + Kd = F,
\end{equation}
where $d$ is the solution vector, $M$, $C$ and $K$ are matrices, and $F$ is a vector, one can make use of methods that damp out spurious high-frequency oscillations such as damped versions of the Newmark method, HHT-$\alpha$, etc.
If one wishes to solve instead a system of the form 
\begin{equation}
M\dot d + Kd = F,
\end{equation}
the obvious choice would be to use a generalized trapezoidal method. However, I am looking for a method that exhibits damping out of spurious high-frequency oscillations, as in the second-order case. I don't need a highly accurate method, but it would preferably be an explicit one.

A:

There is an excellent family of time integration schemes that fit your description called Generalized Single Step Single Solve (GS4).  The original work on the implicit methods for first order systems can be found in [1]. 
Here is the implicit algorithm:
\begin{equation}
C \widetilde{\dot{d}}  + K \widetilde{d}  =  \widetilde{F} 
\end{equation}
where the variables have been approximated between the $n$th and $(n+1)$th timestep as
\begin{align}
& \widetilde{\dot{d}} = \dot{d}_n + \Lambda_6 W_1 (  \dot{d}  _{n+1}-  \dot{d}  _n)
\\
& \widetilde{d}   =  d_n  + \Lambda_4 W_1\Delta t \dot{d}_n + \Lambda_5 W_2\Delta t(  \dot{d}  _{n+1}-  \dot{d}  _n)
\\
& \widetilde{F}  =(1-W_1)q_n+W_1 F_{n+1}
\end{align}
Now one can solve for $ \Delta\dot{d}  =   \dot{d}  _{n+1}-  \dot{d}  _n$ using
\begin{align}
( \Lambda_6 W_1 C &+ \Lambda_5 W_2 \Delta t K)   \Delta \dot{d }
\\
= & -C   \dot{d}  _n - K (   d  _{n} +\Lambda_4 W_1 \Delta t   \dot{d}  _n) 
\\
 &+ (1-W_1)F_n+W_1 F_{n+1}
\end{align}
with the updates
\begin{equation}
\dot{d}  _{n+1} =   \dot{d}  _n +  \Delta \dot{d }  
\end{equation}
\begin{equation}
  d  _{n+1} =   d  _{n} + \lambda_4 \Delta t   \dot{d}  _{n} + \lambda_5 \Delta t  \Delta \dot{d }  
\end{equation}
where
\begin{align}
&\Lambda_4 W_1 = \frac{1}{1+\rho_{\infty}} \nonumber \\
&\Lambda_5 W_2 = \frac{1}{(1+\rho_{\infty})(1+\rho_{\infty}^{s})}\nonumber \\
&\Lambda_6 W_1 = \frac{3+\rho_{\infty}+\rho_{\infty}^{s} - \rho_{\infty}\rho_{\infty}^{s}}{2(1+\rho_{\infty})(1+\rho_{\infty}^{s})}  \\
&W_1 = \frac{1}{1+\rho_{\infty}} \nonumber \\
&\lambda_4 = 1, \quad \nonumber 
\lambda_5 = \frac{1}{1+\rho_{\infty}^{s}} \nonumber
\end{align}
Admittedly, this may seem like a whole load of parameters and nonsense, but all parameters only depend upon two chosen values $(\rho_{\infty}, \rho_{\infty}^{s})$.  The beauty is once you have it programmed you can choose from a whole family of algorithms just by choosing you values for $(\rho_{\infty}, \rho_{\infty}^{s})$.  These parameters correspond to the eigenvalues of amplification matrix of the single DOF problem.  Thus, you can choose the amount of damping (numerical dissipation) simply by choosing them.  Some noteable choices: $(\rho_{\infty}, \rho_{\infty}^{s}) = (1,1)$ gives the Crank-Nicolson method (no damping, not for you) and $(\rho_{\infty}, \rho_{\infty}^{s}) = (0,0)$ gives an algorithm equivalent to the highly dissipative Gear's method (aka two-step backwards difference formula).  Any choice will give you a second-order, unconditionally stable algorithm.
Note the restriction on your choices: $0 \leq \rho_{\infty}^{s} \leq \rho_{\infty} \leq 1$.
Now if you want an explicit algorithm, some algorithms have been developed using the same approach that led to the mess above.  I don't think they have been published any place highly visible yet but early work can be found in a master's thesis here [2].  
The easiest thing to do to get an explicit scheme, with the nice dissipation properties of the above algorithm, is to turn it into a predictor-corrector method.  You, of course, lose the unconditional stability, but you will still have a second-order time integrator.  To do so you can replace the $\widetilde{d}$ above with:
\begin{align}
\widetilde{d}   =  d_n  + \Lambda_4 W_1\Delta t \dot{d}_n
\end{align}
and lump the $C$ matrix, then march to your hearts content.  Everything else stays the same but the restrictions to the $\rho$ parameters above are lifted (they can be anything).  The stability of the algorithm and the amount of dissipation still depends upon this choice.  Hover around 0 and you should be okay.
[1] http://onlinelibrary.wiley.com/doi/10.1002/nme.3228/full
[2] http://conservancy.umn.edu/handle/11299/162393

