Consider two masses 
connected by a single spring 
confined to the 
-axis (
)
![\[L={m\over 2}(\dot{x}_1^2+\dot{x}_2^2)-{m\omega_0^2\over 2}(x_2-x_1-a)^2\]](http://abadonna.physics.wisc.edu/wp-content/ql-cache/quicklatex.com-5037d950c6f6c29a176da29bef23faad_l3.png "Rendered by QuickLaTeX.com")
External forces of dissipation are described by
![\[\mathcal{F}={m\gamma\over 2}(\dot{x}_2-\dot{x}_1)^2\]](http://abadonna.physics.wisc.edu/wp-content/ql-cache/quicklatex.com-b443e4b60860501c7867765aaf327339_l3.png "Rendered by QuickLaTeX.com")
At 
both masses are at stationary equilibrium. Mass one (leftmost) is struck a blow of impulse 
(directed to the right). Find 
for 
in the case 
.
![\[{d\over dt}\Big({\partial L\over \partial \dot{x}_i}\Big)-{\partial L\over \partial x_i}+{\partial \mathcal{F}\over \partial \dot{x}_i}=0, \quad \bm{x}=\left(\begin{array}{c}x_1(t)\\ x_2(t)\end{array}\right)=\left(\begin{array}{c}a\\ b\end{array}\right)e^{i\omega t}=\bm{A}\, e^{i\omega t}\]](http://abadonna.physics.wisc.edu/wp-content/ql-cache/quicklatex.com-dc71f3f5ba6b47943a68f47518ff7126_l3.png "Rendered by QuickLaTeX.com")
![\[\left(\begin{array}{cc} -m\omega^2+k+im\gamma\omega & -k-im\gamma\omega\\ -k-im\gamma\omega & -m\omega^2+k+im\gamma\omega\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right)=0\]](http://abadonna.physics.wisc.edu/wp-content/ql-cache/quicklatex.com-aa5b2399cb0588e47195edf3edbf1027_l3.png "Rendered by QuickLaTeX.com")
![\[(-m\omega^2+k+im\gamma\omega)^2-(-k-im\gamma\omega)^2=0, \quad -m\omega^2+k+im\gamma\omega=\pm(k+im\gamma\omega)\]](http://abadonna.physics.wisc.edu/wp-content/ql-cache/quicklatex.com-4b112519dc14c3ae7d5e3683e6263260_l3.png "Rendered by QuickLaTeX.com")
We can use REDUCE (https://sourceforge.net/projects/reduce-algebra/), one of the most mature open-source CAS projects available (and written by a physicist Tony Hearn) to work out the Lagrange-Euler EOMs and solve the secular equation:
operator x,v;
depend x,t;
depend v,t;
!L:=(m/2)*(v(1)^2+v(2)^2)-(k/2)*(x(2)-x(1))^2;
!F:=m*g*(v(2)-v(1))^2;
EOM1:=df(df(!L,v(1)),t)-df(!L,x(1))+df(!F,v(1));
EOM2:=df(df(!L,v(2)),t)-df(!L,x(2))+df(!F,v(2));
NM1:=sub(x(1)=a*exp(i*w*t), x(2)=b*exp(i*w*t), v(1)=i*w*a*exp(i*w*t), v(2)=i*w*b*exp(i*w*t), EOM1);
NM2:=sub(x(1)=a*exp(i*w*t), x(2)=b*exp(i*w*t), v(1)=i*w*a*exp(i*w*t), v(2)=i*w*b*exp(i*w*t), EOM2);
matrix secular(2,2);
secular(1,1):=df(NM1,a)/exp(i*w*t);
secular(1,2):=df(NM1,b)/exp(i*w*t);
secular(2,1):=df(NM2,a)/exp(i*w*t);
secular(2,2):=df(NM2,b)/exp(i*w*t);
solve({det(secular)},{w});To run REDUCE on the command line, open a terminal and type reduce. Either cut and paste your code into the REPL, or save it as a file (file.red) and instead type reduce -f file.red.
Operator is the generic potentially noncommutative variable that may have dependencies on other variables, which we declare. Clearly df(f,t) is 
, and sub(x=a,f) means substitute 
into the function 
. The output looks like this:
{w=0,
2 2
sqrt( - 2*g *m + k*m)*sqrt(2) + 2*g*i*m
w=------------------------------------------,
m
2 2
- sqrt( - 2*g *m + k*m)*sqrt(2) + 2*g*i*m
w=---------------------------------------------}
m
You can load the “rlfi” package and create LaTeX output. REDUCE is very easy to learn, here is a short collection of all of the basic commands that you will need http://abadonna.physics.wisc.edu/download/reduce-primer/, and REDUCE is very fast (coded in lisp) and runs on the most modest hardware.