Следуя инструкциям по предоставленной вами ссылке, в SymPy это выглядит так:
>>> var('a1 b1 a2 b2 w t')
(a1, b1, a2, b2, w, t)
>>> eq = a1*cos(w*t+b1)+a2*cos(w*t+b2)
>>> cs = cos(w*t), sin(w*t)
>>> xeq=collect(eq.expand(trig=True), cs) # note trig=True
>>> B, A = [xeq.coeff(i) for i in cs]
>>> now=sqrt(A**2+B**2)*cos(w*t-atan2(B,A))