# How to solve 2nd order ODE with matrix parameters in Matlab

1. Dec 3, 2011

### iqjump123

1. The problem statement, all variables and given/known data

I have a frequency equation to solve for the displacement for a spring mass damper truss system, as seen below,

[m]u''+[c]u'+[k]u=f(t),
where m,c,k, are all matrices (2x2), and f(t) is a graph-defined forcing function. I am to use 3 nodes, using the central difference approximation method.

I am to solve for the displacement and the time history of the displacement.

2. Relevant equations

What I have learned so far involved solving matrix equations using [A]{x}={b}

3. The attempt at a solution
I am just confused as to the basic concepts- I should have to be able to decouple this equation, because this is essentially trying to solve a space and time dependent problem.

I am not sure how I can solve for this, and also doing that using the CDM- I would presume that I would somehow solve for the displacement at each time "point", but do i have to use the CDM to solve for the displacements as well? If that is the case, how can I solve the equation that involves matrices and not values?

any help to get me started will be of immense help.

thanks!

2. Dec 4, 2011

### AlephZero

There is something wrong there. If you have 3 nodes your matrices should be 3x3 not 2x2.

Or do you mean "there are 3 nodes in the structure, but one of them is fixed, so it can be eliminated."

Start by understanding how to use CDM for a system with just one variable, for example a mass oscillating on a spring.

You get an equation that says $x_{t+h}$ = something involving $x_t$, $x_{t-h}$, f(t), and k c and m.

For a system with several degrees of freedom, you use exactly the same equation, except x and f are vectors and k c and m are matrices. The only difference is that instead of dividing by a scalar value of m, you multiply by the inverse matrix $m^{-1}$.

CDM is often used when the mass matrix m is diagonal, and inverting a diagonal matrix is easy.

3. Dec 4, 2011

### iqjump123

Hey alephzero,

Thanks so much for your help. It definitely helped in me understanding the approach. I have been consulting some numerical methods handbooks to get myself familiar with when mck are single variables. It seems that I will have to build a tri-diagonal matrix and make it solve for a vector of x values(it seems at the very minimum, I will need x1, x2, xj, xn-2 xn-1 values, with j being the iteration i am solving for, and n-1 being the last iteration. If I am going for more iterations, however, then the xj portion will be lengthened, right?)

However, I am still a bit puzzled in how I can solve this for a matrix of k and c values (since m will be multiplied through by the inverse method)- for example, in row 1 column 1 of the tridiagonal matrix, I will be multiplying by 2+delT*q_1, where q_1 stands for the c matrix corresponding at a certain time. So will my tridiagonal matrix involve another matrix multiplication?

Also, your earlier assumption is correct in that my problem requires a 3 node system. Then the way to solve that is to iterative this code 3 times for the 3 different node and force locations, is that correct?

thanks, as always.

iqjump123

4. Dec 5, 2011

### AlephZero

You seem be getting confused about the the number of variables at different points in space, and the values of their displacement at diffferent times.

If you have a single variable problem like a mass on a spring, the equation of motion equation is mu'' + cu' + ku = f(t)

You convert that into an equation connecting the values of x at three different times, using difference approximations. Call the three times t-h, t, and t+h where h is the time step.

You can approximate u'(t) by [u(t+h) - u(t-h)] / 2h
and u''(t) by [u(t+h) - 2u(t) +u(t-h)] / h^2

So the finite difference version of the equation of motion is

m [u(t+h) - 2u(t) +u(t-h)] / h^2 + c [u(t+h) - u(t-h)] / 2h + k u(t) = f(t)
or
m [u(t+h) - 2u(t) +u(t-h)] + (ch/2) [u(t+h) - u(t-h)] + kh^2 u(t) = h^2 f(t)

As you work through the solution, you already know u(t) and u(t-h) so you can rearrange this to find u(t+h)

[m + ch/2] u(t+h) = (2m - kh^2) u(t) + (-m + ch/2) u(t-h) + h^2 f(t)

If you are starting from t = 0, the initial conditons give you the values of u(0) and u'(0).
To get the calculation started you need the values of u(0) and u(-h). One way to get u(-h) is use a Taylor series approximation and say
u(-h) = u(0) - h u'(0).

So the first time step of the CDM is to calculate u(h) from
[m + ch/2] u(h) = (2m - kh^2) u(0) + (-m + ch/2) u(-h) + h^2 f(0)
For the next time steps, you calculate
[m + ch/2] u(2h) = (2m - kh^2) u(h) + (-m + ch/2) u(0) + h^2 f(h)
[m + ch/2] u(3h) = (2m - kh^2) u(2h) + (-m + ch/2) u(h) + h^2 f(2h)
etc.
You don't need to make one big set of equations for all the time steps.

If you have more than 1 variable (3 in your case), at each step you have to solve a 3x3 set of equations. The equation
[m + ch/2] u(t+h) = (2m - kh^2) u(t) + (-m + ch/2) u(t-h) + h^2 f(t)
is a matrix equation like Au(t+h) = b, where A = m + ch/2. You have to do the matrix multiplcations on the right hand side to calculate the vector b, and then solve the simultaneous equations for u(t+h).

Note, in my first post I said you just needed to invert m. Sorry, I forgot you had the damping matrix c in your model so that wasn't quite right.