Why does this method of solving a system of 2nd order differential equations not work?

[zenterix]

Homework Statement

Consider the system of differential equations

$$\begin{bmatrix} \ddot{x}_1 \\ \ddot{x}_2 \end{bmatrix}=\begin{bmatrix} -\frac{(k_B+k_C)}{m_1} & \frac{k_B}{m_1} \\ \frac{k_B}{m_2} & -\frac{(k_A+k_B)}{m_2} \end{bmatrix}\tag{1}$$

$$\vec{\ddot{x}}=MK^{-1}\vec{x}\tag{2}$$

Relevant Equations

These are the equations of motion for the coupled oscillators in the following picture

I have a question about solving this system.

I (naively, I think) initially did the following

Trial solution: $\vec{x}=e^{\lambda t}\vec{a}$.

Sub this into the system (2)

$$\lambda^2 e^{\lambda t}\vec{a}=M^{-1}Ke^{\lambda t}\vec{a}\tag{3}$$

$$(M^{-1}K-\lambda^2I)\vec{a}=0\tag{4}$$

I then solved this eigenvalue equation for the parameter values

$$m_1=4$$
$$m_2=12$$
$$k_A=78$$
$$k_B=60$$
$$k_C=24$$

and got

$$\lambda_1=-6.37$$

$$\lambda_2=-26.12$$

and associated eigenvectors

$$\vec{a}_1=\begin{bmatrix} 1.025\\1 \end{bmatrix}$$

$$\vec{a}_2=\begin{bmatrix} -2.92 \\ 1 \end{bmatrix}$$

This doesn't seem to be correct in physical terms. These values imply an exponential decay of the positions of the masses towards equilibrium.

The expected result contains oscillations since there is no damping.

Now, it seems to me that the correct path is to transform the second-order system (1) into a first-order system. I know how to do this.

I will do it right after this post (I haven't done it yet because I've done it in another problem and I know it works).

My question is: what is the mathematical reason that the approach depicted above fails?

 

[zenterix]

Here are my thoughts on what the answer to my question might be.

In reviewing what I learned about solving systems of linear differential equations, I noticed that the following results.

1) The matrix exponential $E(t)=e^{tA}$, where $t$ is a real number and $A$ is an $n\times n$ matrix satisfies the differential equation

$$E'(t)=E(t)A=AE(t)$$

2) The only $n\times n$ matrix function satisfying the initial-value problem

$$F'(t)=AF(t),\ \ \ F(0)=B$$

for $-\infty<t<\infty$, is

$$F(t)=e^{tA}B$$

3) $e^{At}=\tilde{\Phi}_0(t)$ where $\tilde{\Phi}_0(t)$ is the normalized fundamental matrix, ie the matrix with linearly independent solutions in each column that satisfy specific initial conditions

$$\tilde{\Phi}_0(0)=I$$

Thus, for the problem in the OP, if we can find two l.i. solutions then we can form this normalized fundamental matrix and we know that we can then form a general solution. General in the sense that every initial value problem solution can be obtained from it.

How to obtain such solutions is actually the starting point of the OP.

One way is from $e^{tA}$ if we can somehow compute it. This is feasible directly only in special cases where a) $e^{tA}$ can be computed in terms of its series definition or b) the matrix $A$ can be diagonalized.

Then there are other specialized methods I have not learned about yet.

But this is irrelevant to the main train of thought I am moving along in this post. The main point is how to obtain two l.i. solutions.

For a first-order system with constant coefficients we do indeed use a trial solution of the form $\vec{x}=e^{\lambda t}\vec{a}$.

I haven't found in my notes a way to solve a system of second-order without transforming it into a linear system by change of variables.

Ok, so in this post I simply went through my current path of thinking without reaching any conclusions. My original question still stands.

 

[TSny]

$$\lambda_1=-6.37$$ $$\lambda_2=-26.12$$

These are the values of $\lambda_1^2$ and $\lambda_2^2$.

 

[kuruman]

I (naively, I think) initially did the following

Also,

• You naively wrote equation (1) to state that a 2×1 matrix on the left is equal to a 2×2 matrix on the right.
• You have three springs, presumably of different spring constants $k_{\!A}$, $k_{\!B}$ and $k_{\!C}.$ Each spring has its own displacement from equilibrium $x_{\!A}$, $x_{\!B}$ and $x_{\!C}$ respectively. How are these three displacements related to the two displacements $x_{1}$ and $x_{2}$ that you show?

It would help if you showed how you obtained your $\mathbf M$ and $\mathbf K$ matrices.

 

[zenterix]

Indeed, last night lying in bed before sleeping I realized that I had calculated the square of $\lambda$.

Continuing on with the calculations, we have

$$\lambda_1^2=-6.37$$

$$\lambda_2^2=-26.12$$

$$\lambda_1=\pm 2.52i$$

$$\lambda_2=\pm 5.11i$$

The solutions $e^{\lambda_1 t}\vec{a}_1$ and $e^{\lambda_2 t}\vec{a}_2$ are complex.

Consider $\lambda_1=\pm 2.52i$. These are complex conjugates. Each of the two eigenvalues gives a complex solution, and each complex solution contains two independent real solutions.

However, the two real solutions of each of these conjugate eigenvalues span the same solution space. Thus, we can choose one of the two conjugate eigenvalues and use the two real solutions it provides to form a general solution for mode $\vec{a}_1$.

The two normal mode solutions are thus

$$\vec{y}_1=\vec{a}_1(c_1\cos{\lambda_1 t}+c_2\sin{\lambda_1 t})$$

$$\vec{y}_2=\vec{a}_2(d_1\cos{\lambda_2 t}+d_2\sin{\lambda_2 t})$$

and the general solution is the sum of these normal mode solutions.

$$\vec{y}=\vec{y}_1+\vec{y}_2$$

The latter is a $2\times 1$ vector representing $x_1$ and $x_2$, the displacements of the two masses in the original problem.

There are four unknown constants $c_1, c_2, d_1,$ and $d_2$ which we can obtain by imposing initial conditions.

 

[zenterix]

Also,

• You naively wrote equation (1) to state that a 2×1 matrix on the left is equal to a 2×2 matrix on the right.
• You have three springs, presumably of different spring constants $k_{\!A}$, $k_{\!B}$ and $k_{\!C}.$ Each spring has its own displacement from equilibrium $x_{\!A}$, $x_{\!B}$ and $x_{\!C}$ respectively. How are these three displacements related to the two displacements $x_{1}$ and $x_{2}$ that you show?

It would help if you showed how you obtained your $\mathbf M$ and $\mathbf K$ matrices.

The first bullet point was not due to naivete. It was just a simple typing omission mistake. Obviously the system is $\vec{\ddot{x}}=M^{-1}K\vec{x}$ as I wrote in (2).

As for the second bullet point, I can show you the calculations I used below.

I applied Newton's 2nd law to each of the masses. Since we are not considering gravity, this involves writing out the spring forces on each mass.

I believe I did this in a perhaps unnecessarily lengthy way, but at this stage it makes me more comfortable to do it this way.

Instead of writing it all out in latex, I'm going to screenshot my handwritten calculations.

Recall the setup

First I write the positions of each of the points $a, b, c, d$, and the wall. I am only considering the vertical direction.

$w$ in these equations represents the width of mass 2. I know that this is unnecessary.

Given these positions, I compute the spring forces.

For example, $F_A=-k_A((x_A-x_{wall})-(x_A-x_{wall})_{eq})$ where $(x_A-x_{wall})_{eq}=-l_A$. The magnitude of the latter is the equilibrium length of spring $A$.

Now that I have the forces, I can use Newton's 2nd law

The OP is based on this last matrix equation.

 

[zenterix]

I would like to show the numerical calculations.

Here is the matrix $M^{-1}K$

$$
\left[\begin{array}{cc}
-\frac{k_{B}+k_{C}}{m_{1}} & \frac{k_{B}}{m_{1}}
\\
\frac{k_{B}}{m_{2}} & -\frac{k_{A}+k_{B}}{m_{2}}
\end{array}\right]
$$

Recall that we want to solve the eigenvalue equation

$$(M^{-1}K-\lambda^2 I)\vec{a}=\vec{0}$$

The $\lambda^2$ eigenvalues are

$$\left[\begin{array}{c}
\frac{-m_{\mathit{1}} k_{A}-m_{\mathit{1}} k_{B}-m_{\mathit{2}} k_{B}-m_{\mathit{2}} k_{C}+\sqrt{k_{A}^{2} m_{\mathit{1}}^{2}+2 k_{A} k_{B} m_{\mathit{1}}^{2}-2 k_{A} k_{B} m_{\mathit{1}} m_{\mathit{2}}-2 k_{A} k_{C} m_{\mathit{1}} m_{\mathit{2}}+k_{B}^{2} m_{\mathit{1}}^{2}+2 k_{B}^{2} m_{\mathit{1}} m_{\mathit{2}}+k_{B}^{2} m_{\mathit{2}}^{2}-2 k_{B} k_{C} m_{\mathit{1}} m_{\mathit{2}}+2 k_{B} k_{C} m_{\mathit{2}}^{2}+k_{C}^{2} m_{\mathit{2}}^{2}}}{2 m_{\mathit{2}} m_{\mathit{1}}}
\\
-\frac{m_{\mathit{1}} k_{A}+m_{\mathit{1}} k_{B}+m_{\mathit{2}} k_{B}+m_{\mathit{2}} k_{C}+\sqrt{k_{A}^{2} m_{\mathit{1}}^{2}+2 k_{A} k_{B} m_{\mathit{1}}^{2}-2 k_{A} k_{B} m_{\mathit{1}} m_{\mathit{2}}-2 k_{A} k_{C} m_{\mathit{1}} m_{\mathit{2}}+k_{B}^{2} m_{\mathit{1}}^{2}+2 k_{B}^{2} m_{\mathit{1}} m_{\mathit{2}}+k_{B}^{2} m_{\mathit{2}}^{2}-2 k_{B} k_{C} m_{\mathit{1}} m_{\mathit{2}}+2 k_{B} k_{C} m_{\mathit{2}}^{2}+k_{C}^{2} m_{\mathit{2}}^{2}}}{2 m_{\mathit{2}} m_{\mathit{1}}}
\end{array}\right]$$

The associated eigenvectors (ie, normal modes) are

$$
\left[\begin{array}{c}
\frac{k_{B}}{\frac{-m_{\mathit{1}} k_{A}-m_{\mathit{1}} k_{B}-m_{\mathit{2}} k_{B}-m_{\mathit{2}} k_{C}+\sqrt{k_{A}^{2} m_{\mathit{1}}^{2}+2 k_{A} k_{B} m_{\mathit{1}}^{2}-2 k_{A} k_{B} m_{\mathit{1}} m_{\mathit{2}}-2 k_{A} k_{C} m_{\mathit{1}} m_{\mathit{2}}+k_{B}^{2} m_{\mathit{1}}^{2}+2 k_{B}^{2} m_{\mathit{1}} m_{\mathit{2}}+k_{B}^{2} m_{\mathit{2}}^{2}-2 k_{B} k_{C} m_{\mathit{1}} m_{\mathit{2}}+2 k_{B} k_{C} m_{\mathit{2}}^{2}+k_{C}^{2} m_{\mathit{2}}^{2}}}{2 m_{\mathit{2}}}+k_{B}+k_{C}}
\\
1
\end{array}\right]
$$

$$
\left[\begin{array}{c}
\frac{k_{B}}{-\frac{m_{\mathit{1}} k_{A}+m_{\mathit{1}} k_{B}+m_{\mathit{2}} k_{B}+m_{\mathit{2}} k_{C}+\sqrt{k_{A}^{2} m_{\mathit{1}}^{2}+2 k_{A} k_{B} m_{\mathit{1}}^{2}-2 k_{A} k_{B} m_{\mathit{1}} m_{\mathit{2}}-2 k_{A} k_{C} m_{\mathit{1}} m_{\mathit{2}}+k_{B}^{2} m_{\mathit{1}}^{2}+2 k_{B}^{2} m_{\mathit{1}} m_{\mathit{2}}+k_{B}^{2} m_{\mathit{2}}^{2}-2 k_{B} k_{C} m_{\mathit{1}} m_{\mathit{2}}+2 k_{B} k_{C} m_{\mathit{2}}^{2}+k_{C}^{2} m_{\mathit{2}}^{2}}}{2 m_{\mathit{2}}}+k_{B}+k_{C}}
\\
1
\end{array}\right]
$$

Now, using the parameter values specified in the OP this all becomes

$$
\left[\left[\begin{array}{c}
- 6.372626865
\\
- 26.12737314
\end{array}\right],
\\
\left[\begin{array}{cc}
1.025474627 & - 2.925474626
\\
1.0 & 1.0
\end{array}\right]\right]
$$

The eigenvalues $\lambda^2$ are the two values on the left matrix, and associated eigenvectors $\vec{a}_1$ and $\vec{a}_2$ are in the columns of the right matrix.

Then

$$\lambda_1=2.52i$$

$$\lambda_2=5.11i$$

where in each case I have chosen the complex conjugate with the positive sign.

We have two complex solutions

$$\vec{y}_{c,1}=e^{\lambda_1 t}\vec{a}_1$$

$$\vec{y}_{c,2}=e^{\lambda_2 t}\vec{a}_2$$

In each of the above complex solutions, we take the real and imaginary parts to form two l.i. real solutions.

A linear combination of these two real solutions is a general solution for one of the normal modes.

The two normal mode solutions are then

$$\vec{y}_1=\vec{a}_1(c_1\cos{\lambda_1 t}+c_2\sin{\lambda_1 t})$$

$$\vec{y}_2=\vec{a}_2(d_1\cos{\lambda_2 t}+d_2\sin{\lambda_1 t})$$

and the general solution is the sum of these normal mode solutions

$$\vec{y}=\vec{y}_1+\vec{y}_2$$

Here are the numerical calculations for these equations

$$
y_{\mathit{1}} = \left[\begin{array}{c}
1.035474627 c_{\mathit{1}} \cos \! \left( 2.524406240 t \right)+
1.035474627 c_{\mathit{2}} \sin \! \left( 2.524406240 t \right)
\\
c_{\mathit{1}} \cos \! \left( 2.524406240 t \right)+c_{\mathit{2}} \sin \! \left( 2.524406240 t \right)
\end{array}\right]
$$
$$
y_{\mathit{2}} = \left[\begin{array}{c}
- 2.925474626 d_{\mathit{1}} \cos \! \left( 5.111494218 t \right)-
2.925474626 d_{\mathit{2}} \sin \! \left( 5.111494218 t \right)
\\
d_{\mathit{1}} \cos \! \left( 5.111494218 t \right)+d_{\mathit{2}} \sin \! \left( 5.111494218 t \right)
\end{array}\right]
$$

$$
y = \left[\begin{array}{c}
1.035474627 c_{\mathit{1}} \cos \! \left( 2.524406240 t \right)+
1.035474627 c_{\mathit{2}} \sin \! \left( 2.524406240 t \right)-
2.925474626 d_{\mathit{1}} \cos \! \left( 5.111494218 t \right)-
2.925474626 d_{\mathit{2}} \sin \! \left( 5.111494218 t \right)
\\
c_{\mathit{1}} \cos \! \left( 2.524406240 t \right)+c_{\mathit{2}} \sin \! \left( 2.524406240 t \right)+
d_{\mathit{1}} \cos \! \left( 5.111494218 t \right)+d_{\mathit{2}} \sin \! \left( 5.111494218 t \right)
\end{array}\right]
$$

At this point, we can impose initial conditions to find the constants.

Here are a few examples of this.

1) Suppose both masses start at equilibrium positions with zero velocity.
Then we get

$$\mathit{constants} = \left\{c_{\mathit{1}} = 0.0,c_{\mathit{2}} = 0.0,d_{\mathit{1}} = 0.0,d_{\mathit{2}} = 0.0\right\}$$

2) Suppose $m_1$ starts at position $-2$ with zero velocity. Then we get
$$\mathit{constants} = \left\{c_{\mathit{1}} = - 0.5049294682,c_{\mathit{2}} = - 0.0,d_{\mathit{1}} = 0.5049294682,d_{\mathit{2}} = 0.0\right\}$$

and the solution $y$ is

$$
\mathit{solution} = \left[\begin{array}{c}
- 0.5228416527 \cos \! \left( 2.524406240 t \right)- 1.477158347 \cos \! \left( 5.111494218 t \right)
\\
- 0.5049294682 \cos \! \left( 2.524406240 t \right)+
0.5049294682 \cos \! \left( 5.111494218 t \right)
\end{array}\right]
$$

and if we plot this solution (again, the first row represents displacement of mass 1 and the second row the displacement of mass 2) then we get

 

[zenterix]

Here is another scenario: both $m_1$ and $m_2$ start at position $-2$ with zero velocity.

$$\mathit{constants} = \left\{c_{\mathit{1}} = - 1.982087815,c_{\mathit{2}} = - 0.0,d_{\mathit{1}} = - 0.01791218455,d_{\mathit{2}} = - 0.0\right\}$$

$$
\mathit{solution} = \left[\begin{array}{c}
- 2.052401641 \cos \! \left( 2.524406240 t \right)+

0.05240164140 \cos \! \left( 5.111494218 t \right)
\\
- 1.982087815 \cos \! \left( 2.524406240 t \right)-

0.01791218455 \cos \! \left( 5.111494218 t \right)
\end{array}\right]$$

The red curve represents mass 1.

Not sure, intuitively, why mass 1 has a larger amplitude.

If we make the masses the same at $4$ then we get

This looks weird, but if I try to simulate something similar using MIT's coupled oscillators applet I do get something somewhat similar

Note that the setup is not exactly the same in the MIT applet, but it is similar.

With my calculations and similar parameters to the above I get

Finally, let's check what happens when the masses are the same and the spring constants are all the same.

 

[kuruman]

This

assumes that the bar stays parallel to itself at all times. The setup as you described it is not symmetric about the vertical that passes through the midpoint of the horizontal bar. You have $m_2$ on the left and nothing on the right. Also the springs have all different spring constants.

This means that the instantaneous forces at $c$ and $d$ will be in general different generating a net torque about the middle of the rod. As a result, there will be a frequency of torsional oscillations of the bar to consider that will be coupled to its linear modes of oscillations of $m_1$ and $m_2$.

To see what's involved, you may wish to try solving a simpler problem. Get rid of $m_2$ and attach each end of the bar (moment of inertia $I=\frac{1}{12}ML^2$ about its CM) to springs of the same equilibrium length but unequal constants $k_A$ and $k_B$. Displace each end by different amounts from the equilibrium position and release from rest. The problem that you posted is more difficult than this because of $m_2.$

 

[zenterix]

@kuruman

That is true. Please see the original problem here. It is the last problem in the problem set. I wonder what they have in mind regarding the issue of torque.

There is the assumption that "the unstretched lengths of spring A and B equal the unstretched length C", but this doesn't seem to solve the issue of torque on mass 1.

There is also the assumption that "the masses are constrained to move only in the vertical direction". If mass 1 is somehow impeded from rotating and can only move up or down then I think maybe we can say the calculations I did are correct?

I actually just realized that they are including gravity in the problem. I think this only changes the equilibrium positions. I will add gravity.

 

[kuruman]

It looks like the author of the problem did not consider the possibility of torque. Under the assumption that the torsional motion is ignored, your frequencies look correct. I checked that they reduce to the correct expression s when $m_2\rightarrow 0.$ As you say, adding gravity will not change the frequencies. Just redefine your $l_a,~l_b,~l_c$ to be the lengths of the springs at equilibrium.

 

[zenterix]

Here is what I think happens when we add gravity to the problem.

The equations of motion become

$$F_a+F_b-m_2g=m_2\ddot{x}_2$$

$$F_c+F_d-m_1g=m_1\ddot{x}_1$$

which become

$$\ddot{x}_2=\frac{k_B}{m_2}x_1-\frac{(k_A+k_B)}{m_2}x_2-g$$

$$\ddot{x}_1=-\frac{(k_B+k_C)}{m_1}x_1+\frac{k_B}{m_1}x_2-g$$

The previous calculations showed how to obtain the solution to the homogeneous system, which does not have the $-g$ terms in the equations.

For a first-order linear system we can use variation of parameters to obtain a solution to the inhomogeneous equation.

What is the method that is commonly used to solve for this solution when we have a second-order system?

 

[zenterix]

This post is just to illustrate the problem further.

The following are plots of the two normal mode solutions for the same initial conditions $x_1(0)=2$ and $x_1'(0)=0$. That is, these are the solutions in which the first mass starts at position 2 with zero velocity.

The normal mode solutions arise from the matrix $M^{-1}K$ and the fact that every solution is an exponential.

$M^{-1}K$ represents the specifics of the physical system.

The square root of an eigenvalue appears as a factor in the exponent of an exponential solution and represents an angular frequency.

If such a square root is complex then we have oscillations.

The associated eigenvectors are called normal modes and they affect the amplitudes when we have oscillations.

In the case of the system in the OP, one of the normal modes has components with the same sign and the other normal mode has components with opposite sign.

Because each normal mode is associated with a single angular frequency, a solution involving just the first normal mode has the masses moving with same angular frequency and phase and in the same direction at all times.

A solution involving just the second mode has the masses moving with the same angular frequency and phase (albeit different from the frequency of the first mode) but in opposite directions.

One question I was asking myself just now is how to reconcile the fact that we have two different solutions for the same initial conditions (given what was written in post #2 above).

First of all, the system we are solving is second-order. The theorems in post #2 are for a first-order system.

If we transform the 2nd-order system into first-order it becomes a $4\times 4$ system.

For that system, there is a unique solution $e^{tA}$ where

$$A=\begin{bmatrix} 0_2&I_2\\ M^{-1}K& 0_2 \end{bmatrix}$$

is a $4\times 4$ matrix.

$e^{tA}$ equals the normalized fundamental matrix of this system, which means it four independent solutions in its columns.

We can find these solutions by calculating eigenvalues and eigenvectors. There are four eigenvalues for $A$ and they are composed of two pairs of complex conjugates.

Since these eigenvalues are exponents in an exponential, they represent angular frequency when we obtain real solutions. Solutions obtained from a pair of conjugates represent the same set of solutions.

Thus, we take one complex solution from each conjugate pair and form two independent real solutions with the same angular frequency.

Since we do this twice, we obtain four total real solutions.

Linear combinations of the solutions with the same angular frequency represent a normal mode.

 

[kuruman]

What is the method that is commonly used to solve for this solution when we have a second-order system?

It's simple substitution in this case. Consider a simple vertical mass-spring system with gravity.
First you write $~\ddot x=-\dfrac{k}{m}x-g$

Then you let $\xi=x-\dfrac{mg}{k} \implies \ddot {\xi}=\ddot x$ and $~x=\xi+\dfrac{mg}{k}$

so that the equation of motion to solve becomes
$\ddot {\xi}=-\dfrac{k}{m}\left(\xi+\dfrac{mg}{k}\right)-g\implies \ddot {\xi}=-\dfrac{k}{m}\xi.$

Same equation to solve except that variable $\xi$ is the displacement from the equilibrium position which is at $x_{eq}=\dfrac{mg}{k}.$