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

[iqjump123]

Homework Statement

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.

Homework Equations

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

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!

 

[AlephZero]

[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.

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."

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.

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.

 

[iqjump123]

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.

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

 

[AlephZero]

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?)

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.