System of linear second order differential equations

[matematikawan]

I think I know how to solve
\frac{d\vec{x}}{dt}= A \vec{x}

where A is a constant nXn matrix. We just compute the eigenvalues and the corresponding eigenvectors.

But how do we solve
\frac{d^2\vec{x}}{dt^2}= A \vec{x}

Can we say straight away that the solution is (following that of one dependent variable)
\vec{x}(t) = exp(-Mt) \vec{x}_1+ exp(Mt) \vec{x}_2

where M²=A and x₁ and x₂ are constant vectors.

 

[tiny-tim]

Hi matematikawan!

Try putting y = dx/dt + Mx

 

[matematikawan]

This is great. How did you came out with the substitution y = dx/dt + Mx tiny-tim ?. I have try searching the internet for system of DE. Most of them are about first order DE.
With your substitution I obtained
\frac{d\vec{y}}{dt}=M\vec{y}

Using http://en.wikipedia.org/wiki/Matrix_exponential"
\vec{y} = exp(Mt) \vec{c_1}

\frac{d\vec{x}}{dt} = -M\vec{x} + exp(Mt) \vec{c_1}

\vec{x}=exp(-Mt)\vec{c_2} + exp(-Mt)\int exp(2Mt) \vec{c_1}

\vec{x}=exp(-Mt)\vec{c_2} + exp(-Mt) (2M)^{-1} exp(2Mt) \vec{c_1}


So not exactly

\vec{x}(t) = exp(-Mt) \vec{x}_1+ exp(Mt) \vec{x}_2

as I predicted.

 

[matematikawan]

Or it could possible that

\vec{x}=exp(-Mt)\vec{c_2} + exp(-Mt) exp(2Mt) (2M)^{-1}\vec{c_1}

Then it is of the form

\vec{x}(t) = exp(-Mt) \vec{x}_1+ exp(Mt) \vec{x}_2

 

[tiny-tim]

Hi matematikawan!

(btw, we usually write the constant before the exponential, and we usually use A and B if we're expecting more than one )

This is great. How did you came out with the substitution y = dx/dt + Mx tiny-tim ?. I have try searching the internet for system of DE. Most of them are about first order DE.
With your substitution I obtained
\frac{d\vec{y}}{dt}=M\vec{y}

That's right …

now write it dx/dt + Mx = Be^Mt,

and use the particular solution and general solution method (do you know what that is?) …

the g s gives you Ce^-Mt, and the p s of course is (B/M - 1)e^Mt, or just Be^Mt.

And I got it from the "D" method …

we write D instead of d/dt, so the original equation is (D² - A)x = 0,

then (D + A)(D - A)x = 0,

so either (D + A)x = 0 or (D - A)x = 0 …

this works for any polynomial in D (it's called the characteristic polynomial of the original equation) …

google "characteristic polynomial", or see the recurrence relation page in the PF Library for some details. ​

 

[matematikawan]

That is a clever trick. I have seen in some old DE textbooks where they use similar 'D operator' technique to solve ODE. But that works only for DE with one dependent variable.

In this case of system of DE

\frac{d^2\vec{x}}{dt^2}= A \vec{x}

I'm quite lucky to be able to find M such that M^2=A. To find M in MATLAB we just use the command sqrtm(A).

How about these system of DE
\frac{d^3\vec{x}}{dt^3}= A \vec{x}

or
\frac{d^2\vec{x}}{dt^2} + \frac{d\vec{x}}{dt}= A \vec{x}

How are we going to 'factorise' the so called characteristic polynomial ?

 

[tiny-tim]

How are we going to 'factorise' the so called characteristic polynomial ?

(erm … it's not "so-called" … that really is its name!)

For the first one, use the cube root of A (I don't know whether there's a MATLAB command for that: if not, is there one for lnA or expA?).

The second is easy: the roots are -(1/2)I ±√((1/4)I + A²).

 

[matematikawan]

Oh I'm sorry. :shy: I didn't mean to offend anyone. May be I'm not using proper words there. English is not my first language. So it is characteristic polynomial. Really sorry.

But by the way when I Google 'characteristic polynomial'. This term gives the meaning as the polynomial whose roots are the eigenvalues of a given matrix A, i.e. det(A-aI)=0.

In MATLAB we use the command expm(A) for matrix exponential. I'm not aware of commands for cube root or log of A. I'm trying to solve second order pde numerically. I discretise the variables x and t. Solution for the third order system of DE is not that important for me at the moment compare to that of second order.
It looks like that we can factorise the characteristic quadratic by completing the square.

D^2+D-A =\left(D + \frac{1}{2}I\right)^2 - \frac{1}{4} I - A
= \left(D + \frac{1}{2}I - \sqrt{\frac{1}{4}I+A} \right) \left(D + \frac{1}{2}I + \sqrt{\frac{1}{4}I+A} \right)

Thank you for your help.

 

[D H]

Hi matematikawan!

Try putting y = dx/dt + Mx

Alternatively, form a state vector comprising \vec x and d\vec x/dt.

<br /> \vec u = \left [\vec x \\ d\vec x/dt \right]<br />

Then

<br /> \dot{\vec u} =<br /> \bmatrix \boldsymbol 0_{N\times N} &amp; \boldsymbol 1_{N\times N} \\<br /> \boldsymbol A &amp; \boldsymbol 0_{N\times N} \endbmatrix \vec u<br />

where \boldsymbol 1 is the NxN identity matrix and \boldsymbol 0 is the NxN zero matrix. Denoting

\boldsymbol B \equiv \bmatrix<br /> \boldsymbol 0_{N\times N} &amp; \boldsymbol 1_{N\times N} \\<br /> \boldsymbol A &amp; \boldsymbol 0_{N\times N} \endbmatrix

This is of course a constant matrix. Then

\dot{\vec u} = \boldsymbol B \vec u

and tada! A linear first order equation with a constant matrix.

 

[matematikawan]

Alternatively, form a state vector comprising \vec x and d\vec x/dt.

<br /> \vec u = \left [\vec x \\ d\vec x/dt \right]<br />

Then

<br /> \dot{\vec u} =<br /> \bmatrix \boldsymbol 0_{N\times N} &amp; \boldsymbol 1_{N\times N} \\<br /> \boldsymbol A &amp; \boldsymbol 0_{N\times N} \endbmatrix \vec u<br />

Of course, we could then solve the equation with Runge-Kutta method, thus avoiding the need for matrix exponential. Although we have double the matrix size but their entries are sparse. Anyway I'm looking for numerical answer. Analytical answer may be suitable for theoretical purposes.

Thanks D H, especially for formulating the equation with nice block matrices. Is there any chance of having B as
\boldsymbol B \equiv \bmatrix<br /> \boldsymbol A &amp; \boldsymbol 0_{N\times N} \\<br /> \boldsymbol 0_{N\times N} &amp; \boldsymbol 1_{N\times N} \endbmatrix

I have try a few permutation but without success.