algebraic manipulation - How to linearize an expression automatically?

I would like to automatically linearize some long equations in the scope of variational calculus. Here follows an example of what I need to do :

Given two variables $a_1 = q_1 + \delta q_1$ and $a_2 = q_2 + \delta q_2$ and a product

$${a_1}^2\, a_2 = {q_1}^2 q_2 + 2q_1q_2\delta q_1 + q_2{\delta q_1}^2 + {q_1}^2\delta q_2 + 2q_1\delta q_1 \delta q_2 + {\delta q_1}^2 \delta q_2$$

I would like to eliminate any variable preceded by the $\delta$ symbol which power is superior to 1 (make it equal to zero), and any product of two variables preceded by the $\delta$ symbol (make the product equal to zero also). So as to obtain :

$${a_1}^2\, a_2 = {q_1}^2 q_2 + 2q_1q_2\delta q_1 + {q_1}^2\delta q_2$$

I first tried the Assumptions options while expanding :

a1 = Subscript[q, 1] + Subscript[\[Delta]q, 1];
a2 = Subscript[q, 2] + Subscript[\[Delta]q, 2];
Expand[a1^2*a2, Assumptions -> Subscript[\[Delta]q, 1]^2 = 0]

Which returned the following :

Set::write: Tag Rule in Assumptions->Subsuperscript[\[Delta]q, 1, 2] is Protected. >>
(Subscript[q, 1] + Subscript[\[Delta]q, 1])^2 (Subscript[q, 
   2] + Subscript[\[Delta]q, 2])

Of course it didn't work. Truth is that I don't know how to start this... Does someone has any ideas?

I also tried :

a1 = Subscript[q, 1] + Subscript[\[Delta]q, 1];
a2 = Subscript[q, 2] + Subscript[\[Delta]q, 2];
b = Expand[a1^2*a2];
Assuming[Subscript[\[Delta]q, 1]^2 == 0, b]

which didn't work either and returned :

\!\(

\*SubsuperscriptBox[\(q\), \(1\), \(2\)]\ 
\*SubscriptBox[\(q\), \(2\)]\) + 
 2 Subscript[q, 1] Subscript[q, 2] Subscript[\[Delta]q, 1] + 
 Subscript[q, 2] 
\!\(\*SubsuperscriptBox[\(\[Delta]q\), \(1\), \(2\)]\) + \!\(
\*SubsuperscriptBox[\(q\), \(1\), \(2\)]\ 
\*SubscriptBox[\(\[Delta]q\), \(2\)]\) + 
 2 Subscript[q, 1] Subscript[\[Delta]q, 1] Subscript[\[Delta]q, 
  2] + \!\(
\*SubsuperscriptBox[\(\[Delta]q\), \(1\), \(2\)]\ 

\*SubscriptBox[\(\[Delta]q\), \(2\)]\)

Answer

I'd use Series :

f[a1_, a2_] = a1^2 a2;

(Series[f[q1 + \[Epsilon] dq1, q2 + \[Epsilon] dq2], {\[Epsilon], 0, 1}] // Normal) 
  /. \[Epsilon] -> 1
(* dq2 q1^2 + 2 dq1 q1 q2 + q1^2 q2 *)