Calculating uncertainties for a final result Say you are dividing 2 times with uncertainties: 
$$\frac{t_1}{t_2} ~=~ \frac{0.551s \pm 0.002s}{ 0.712s \pm 0.002s}.$$
After doing the calculations you get:
$$\frac{t_1}{t_2} ~=~ 0.774 \pm  \ldots?$$   
How do you calculate your uncertainty with the two uncertainties $\pm0.002s$?
 A: The ultimate answer is the JCGM 100:2008 guide followed by most of the metrology institutes around the world. The specific chapter on combining uncertainties is Chapter 5. 
Specifically, for a two-variable function $f(t_1, t_2)$ of two random variables, Eq. (16) of Section 5.2.2. gives $$\Delta f^2= \left (\frac{\partial f}{\partial t_1} \right )^2 \Delta t_1^2+
\left (\frac{\partial f}{\partial t_2} \right )^2 \Delta t_2^2+ 2 \frac{\partial f}{\partial t_1} \frac{\partial f}{\partial t_2} \Delta t_1 \Delta t_2 r(t_1,t_2)$$
where $-1<r(t_1,t_2)<+1$ is the correlation coefficient between $t_1$ and $t_2$. 
Even more specifically, for $f(t_1,t_2)=t_1/t_2$ 
this procedure gives
$$\frac{\Delta f^2}{f^2} = \frac{\Delta t_1^2}{t_1^2}+ \frac{\Delta t_1^2}{t_2^2}-2 r(t_1,t_2) \frac{\Delta t_1}{t_1} \frac{\Delta t_2}{t_2} $$
If $t_1$ and $t_2$ are uncorrelated, we recover the usual formula. In the other extreme, when $t_1=t$ and $t_2=t$ (both are identical to a single random variable) we
get
$$\Delta \frac{t_1}{t_2} = \Delta 1 = \Delta t_1^2+\Delta t_2^2 - 2\Delta t_1 \Delta t_2=(\Delta t_1-\Delta t_2)^2=0$$ 
as we should.
A: Here's the general derivation of the commonly used, and often (but not always) valid, uncertainty propagation formula for independent small Gaussian errors.
$\newcommand{\bbv}[1]{\mathbf{#1}}$
Consider a quantity $y$, calculated from measured quantities $\bbv{x}$
$$
  y + \Delta{y} = f(\bbv{x}+\Delta\bbv{x}) = f(x_1+\Delta{x}_1,x_2+\Delta{x}_2,\ldots,x_n+\Delta{x}_n)
$$
Note that I don't write $\pm$ symbols: it's the actual, positive or negative, deviation of a single measured value from the true values that I mean. The $\Delta x$ are not the measurement uncertainties, but randomly-distributed physical quantities that have somehow found their way into our measurement.
Assume now the deviations to be


*

*independent, i.e. stemming from measurements not directly coupled

*small


where by small I mean, we can Taylor expand the function:
$$
  y + \Delta{y} \approx f(\bbv{x}) + \Delta\bbv{x}\cdot\bbv\nabla{f(\bbv{x})}
  = f(\bbv{x}) + \Delta{x_1}\cdot\frac{\partial}{\partial{x_1}}{f} + \Delta{x_2}\cdot\frac{\partial}{\partial{x_2}}{f} + \ldots
$$
with only the terms I've writte there. Since $f(\bbv{x})=y$, it's just
$$
  \Delta{y} = \Delta{x_1}\cdot\partial_1{f} + \Delta{x_2}\cdot\partial_2{f} + \ldots
$$
(that's not the result yet!)
That the errors are independent means that there exists a probability density function for each of the measured quantities $\psi_i(\Delta{x_i})$. What we're interested in is the corresponding PDF of the result uncertainties
$\newcommand{\ttd}{\mathrm{d}}$
$$\begin{aligned}
  \phi(\Delta{y}) =& \text{probability density of }\bigl(f(\bbv{x}+\Delta\bbv{x})=y+\Delta{y}\bigr)
  \\=& \int_{\{(\Delta x_1,\ldots,\Delta x_n)\in\mathbb{R}^n:f(\bbv{x}+\Delta\bbv{x})=y+\Delta{y}\}}\!\!\!\!\!\!\!\!\!\!\!\!\ttd\!\Delta\!x_1\!\cdot\!\psi_1(\Delta x_1)\ \ttd\!\Delta\!x_2\!\cdot\!\psi_2(\Delta x_2)\ \ldots
\end{aligned}$$
integrating over an $(n-1)$-dimensional manifold (this strictly needs to be written in a less ambiguous way, but it's ok for our purpose).
Now put in the Taylor series approximation
$$\begin{aligned}
  \phi(\Delta{y}) =& \int_{\{(\Delta x_1,\ldots,\Delta x_n):\Delta{y} = \Delta{x_1}\cdot\partial_1{f} + \Delta{x_2}\cdot\partial_2{f}+\ldots\}}\!\!\!\!\!\!\!\!\!\!\!\!\ttd\!\Delta\!x_1\!\cdot\!\psi_1(\Delta x_1)\ \ttd\!\Delta\!x_2\!\cdot\!\psi_2(\Delta x_2)\ \ldots
\end{aligned}$$
and note that this is just the convolution
$$
  (f_1\star f_2\star\ldots\star f_n)(y) = \int_{\{(x_1,x_2,\ldots,x_n)\colon\sum_{x_n}=y\}}\!\!\!\!\!\!\!\!\!\!\!\!\ttd x_1\!\cdot\!f_1(x_1)\ \ttd x_2\!\cdot\!f_2(x_2)\ \ldots
$$
of the functions $\psi_i$. If these are well-behaved in a very weak sense, we can apply the convolution theorem
$$
  \bigl(\mathrm{FT}(f_1\star f_2\star\ldots\star f_n)\bigr)(k) = \prod_i\mathrm{FT}(f_i)(k)
$$
(depending on your definition of the Fourier transform, with some scaling factor). In our case,
$$
  \mathrm{FT}(\phi)(k) = \prod_i\mathrm{FT}\left(\lambda x.\psi_i\Bigl(\frac{x}{\partial_if}\Bigr)\right)(k)
$$
To actually compute this now, we need to assume some specific shape of the functions $\psi_i$, i.e. of the "shape" of the deviations in our measurements. And for such things, a normal distribution is very often a good bet:
$$
  \psi_i(\Delta{x}) = N_i\cdot \exp\Bigl(\frac{-\Delta\!x^2}{2\cdot\sigma{x}_i{}^2}\Bigr).
$$
(Why did nobody correct me? I'd written complete rubbish here!)
We're not really interested in the normalization constant $N_i$, but the standard deviation $\sigma{x}_i$ is highly relevant: this is the actual uncertainty of the measured quantity $x_i$, the thing that's normally meant when writing $\Delta{x_i}$.
The Fourier transform of such a Gaussian bell curve is, conveniently, again a normal distribution:
$$
  \mathrm{FT}(\psi_i)(k) \propto \exp\Bigl(\frac{-k^2\cdot \sigma{x}_i{}^2}{2}\Bigr).
$$
We actually need "stretched" versions of the functions,
$$
  \mathrm{FT}\bigl(\lambda x.\psi_i(\tfrac{x}{s})\bigr)(k)
  = \mathrm{FT}\left(\exp\Bigl(\frac{-\frac{\Delta\!x^2}{s^2}}{2\cdot\sigma{x}_i{}^2}\Bigr)\right)(k)
  \propto e^{\frac{-k^2\cdot s^2\cdot\sigma{x}_i{}^2}{2}}.
$$
Using this, we obtain
$$\begin{aligned}
  \mathrm{FT}(\phi)(k) \propto \prod_i\exp\Bigl(\frac{-k^2\cdot (\partial_if)^2\cdot \sigma{x}_i{}^2}{2}\Bigr)
  \\=& \exp\Bigl(\frac{-k^2}2\sum_i(\partial_if\cdot \sigma{x}_i)^2\Bigr)
\end{aligned}$$
which is yet again a Gaussian function, so we can easily back-transform it,
$$\begin{aligned}
  \phi(\Delta{y}) \propto \mathrm{IFT}\left(\exp\Bigl(\frac{-k^2}2\sum_i(\partial_if\cdot \sigma{x}_i)^2\Bigr)\right)(\Delta{y})
  \\\propto& \exp\Bigl(\frac{-\Delta\!y^2}{2\sum_i(\partial_if\cdot \sigma{x}_i)^2}\Bigr).
\end{aligned}$$
Here, we can associate $\sum_i(\partial_if\cdot \sigma{x}_i)^2$ with the (squared) standard deviation of this normal distribution. So, summarizing, the uncertainty of the quantity
$$
  y = f(\bbv{x})
$$
is
$$
  \sigma{y} = \sqrt{\sum_i\Bigl(\frac{\partial f}{\partial x_i}\cdot \sigma{x_i}\Bigr)^2}
$$
and that's it!
A: $\newcommand{\u}[0]{\pm\Delta} \newcommand{n}[1]{(#1\u #1)} \newcommand{\c}[1]{\frac{\Delta #1}{#1}}$
I'm using the notation $\u a$ to denote uncertainity of a variable $a$
Good method (worst-case error)
This method brings out the error assuming worst-case scenarios.
Let's say you have a bunch of variables $a_i$, and another bunch $b_i$. Lets multiply all the $a$s and divide by the $b$s. So we have something of the sort $$c=\frac{\n{a_1}\n{a_2}\cdots\n{a_i}\cdots\n{a_n}}{\n{b_1}\n{b_2}\cdots\n{b_i}\cdots\n{b_m}}$$
Now, the procedure to calculate uncertainity is simple: First maximise the error in the positive direction by taking $+$ in all numerator terms($\u a_i\leftarrow +\Delta a_i$), and $-$ in the denominator terms ($\u b_i\leftarrow +\Delta b_i$). Calculate the value as $c_+$. Now, do the opposite to maximise the error in the negative direction by swapping the signs, calculate as $c_-$. Now, take the average of the two $c=\frac{c_+ + c_-}{2}$, and $\Delta c= c_+ -c$.
Wrong method
This must not be confused with the following method for approximately computing perturbations. This method gives wrong error bars when applied (with signs of the Deltas) deleted. It,s OK to use for approximations.
First calculate mean value of $c$ by removing all the uncertainity terms and just dividing. Now, taking positive values of all uncertainities, use:
$$\c{c}=\c{a_1}+\c{a_2}+\cdots+ \c{a_i}+\cdots+\c{a_n}+\c{b_1}+\c{b_2}+\cdots+ \c{b_i}+\cdots+\c{b_n}$$
One can see here that if there is a square term $p^2$, the corresponding term in the above formula can be written as $2\c{p}$. This generalised to fractional exponents &c. Generally, for a $p^r$ term, the corresponding term in the above formula is $|r|\c{p}$.
If you're interested, that formula can be easily derived by taking the logarithm of the product, splitting up into a sum, differentiating, and replacing $\frac{\mathrm{d}a_i}{a_i}$ with $\c{a_i}$ (thus the small-uncertainities only)
If the two variables have a relation between them ($t_1=2t_2$ or $t_1=t_2+c$ &c), then use @Slavik's method posted below. For a larger number of variables, the method becomes icky unless you do it with matrices and Taylor expansions.
