Thread Subject: Simple ODE: Linearly Damped Spring Mass System

In Matlab, How can I solve a Linear Second Order Homogenous ODE in Matrix form, i.e.

M x’’ + C x’ + k x = 0
Where M, C, K are 6x6 complex matrices, and x is a function of time.
and some initial conditions, i.e. x(0)=1.

I cannot find references that address this problem in the matrix form. Any tips will be appricated. Thanks,

"Hydroman S" <amirgsalem@gmail.com> wrote in message <gkjq5e$ph5$1@fred.mathworks.com>...
> In Matlab, How can I solve a Linear Second Order Homogenous ODE in Matrix form, i.e.
>
> M x’’ + C x’ + k x = 0
> Where M, C, K are 6x6 complex matrices, and x is a function of time.
> and some initial conditions, i.e. x(0)=1.
>
> I cannot find references that address this problem in the matrix form. Any tips will be appricated. Thanks,

Why do you need any reference? take a look of doc ODE?? (especially Mass matrix part), your system can be solved by those functions without any further manipulation.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gkjueb$1ts$1@fred.mathworks.com>...
> "Hydroman S" <amirgsalem@gmail.com> wrote in message <gkjq5e$ph5$1@fred.mathworks.com>...
> > In Matlab, How can I solve a Linear Second Order Homogenous ODE in Matrix form, i.e.
> >
> > M x’’ + C x’ + k x = 0
> > Where M, C, K are 6x6 complex matrices, and x is a function of time.
> > and some initial conditions, i.e. x(0)=1.
> >
> > I cannot find references that address this problem in the matrix form. Any tips will be appricated. Thanks,
>
> Why do you need any reference? take a look of doc ODE?? (especially Mass matrix part), your system can be solved by those functions without any further manipulation.
>
> Bruno

One more remark, I guess your variable x is (6 x 1) vector.

Bruno

Thanks so much Bruno, and sorry for the late reply.

so here is my attempt in solving the problem using the eig. value problem:


m = rand(6,6); c = rand(6,6); k = rand(6,6); % note these are not the actual values, the are to illustrate my problem here only

w=7; tlim=[0,100]; nt=400;

n=size(m,1); t=linspace(tlim(1),tlim(2),nt);
% Initial conditions
y0=[1;1;1;1;1;1]; v0=zeros(6,1); y0=0*y0;

% eigenvalues and eigenvectors for the homogeneous solution
A=[zeros(n,n), eye(n,n); -m\[k, c]];
[U,lam]=eig(A); [lam,j]=sort(diag(lam)); U=U(:,j);

% homogeneous solution
U=U*diag(U\[y0; v0]);
yh=real(U(1:n,:)*exp(lam*t));

Problems:

when I use the real part of complex m, c & k, the eig. values are:
lam =

        0
        0
        0
  -0.0023
  -0.1068
  -0.1068
  -0.1418 - 0.7198i
  -0.1418 + 0.7198i
  -0.1418 - 0.7198i
  -0.1418 + 0.7198i
  -0.0973 - 0.8748i
  -0.0973 + 0.8748i

Problem: the above method does not like zero eig. values. so my second option is to use ode45 and solve the problem numerically, but I'm not sure how to find the homogeneous solution using this numerical function? any tips will be great.

> Problem: the above method does not like zero eig. values. so my second option is to use ode45 and solve the problem numerically, but I'm not sure how to find the homogeneous solution using this numerical function? any tips will be great.
===============================================================================================================================================================


One more thing: In one reference I found on the file exchange section of this site which deals with a similar problem, but with a forcing function F(t) such that the equation of motion will have a homogenous + a particular solution.

M*X"+C*X'+K*X=F(t)

 A typical call to strdynrk function is:
 m=eye(3,3); k=[2,-1,0;-1,2,-1;0,-1,2];
 c=.05*k; x0=zeros(3,1); v0=zeros(3,1);
 t=linspace(0,10,101);
 [t,x,v]=strdynrk(t,x0,v0,m,c,k,'func');

global Mi C K F n n1 n2
Mi=inv(m); C=c; K=k; F=functim;
n=size(m,1); n1=1:n; n2=n+1:2*n;
[t,z]=ode45(@sde,t,[x0(:);v0(:)]);
x=z(:,n1); v=z(:,n2);
 
%================================

function zp=sde(t,z)
% zp=sde(t,z)
global Mi C K F n n1 n2
zp=[z(n2); Mi*(feval(F,t)-C*z(n2)-K*z(n1))];

%================================

function f=func(t)
% f=func(t)
% This is an example forcing function for
% function strdynrk in the case of three
% degrees of freedom.
f=[-1;0;2]*sin(1.413*t);

My question then becomes, how to get the homogenous solution zh instead of zp, simply making f = [0;0;0] will not do the job, right?

Sorry to bring this up again, but I was hoping that maybe anyone who is familiar with solving ODE’s numerically in Matlab could help me finalize the problem of finding the homogenous solution using ODE45 based on the equation below in the matrix form:

M*X"+C*X'+K*X=0

> Sorry to bring this up again, but I was hoping that
> maybe anyone who is familiar with solving ODE’s
> numerically in Matlab could help me finalize the
> problem of finding the homogenous solution using
> ODE45 based on the equation below in the matrix form:
>
> M*X"+C*X'+K*X=0
>

Your system can be written as a first-order system
as
X' = Y
M*Y' = -K*X-C*Y
or (if M is regular)
X' = Y
Y' = -M^(-1)*K*X-M^(-1)*C*Y
Write this system as
Z' = A*Z
with Z=(X,Y) and A = [0,I ; -M^(-1)*K, -M^(-1)*C]
Then (X,Y)=exp(A*t)*(X_0,Y_0) is the solution of
the homogenous equation with exp being the matrix
exponential.

Take a look at
http://www.mathworks.com/access/helpdesk/help/toolbox/symbolic/index.html?/access/helpdesk/help/toolbox/symbolic/expm.html&http://www.google.de/search?hl=de&q=MATLAB+%26+symbolic+matrix+exponential&meta=
for the calculation of the matrix exponential.

Best wishes
Torsten.

"Hydroman S" <amirgsalem@gmail.com> wrote in message
news:glpnvt$pa8$1@fred.mathworks.com...
> Sorry to bring this up again, but I was hoping that maybe anyone who is
> familiar with solving ODE’s numerically in Matlab could help me
> finalize the problem of finding the homogenous solution using ODE45 based
> on the equation below in the matrix form:
>
> M*X"+C*X'+K*X=0
>

Convert the 2nd order ODE into a system of 1st order ODEs [massMatrix*y' =
f(t, y)] and use ODE45 to solve this system. If you're not sure how to
generate the system of ODEs from the single 2nd order ODE, go to the support
website:

http://www.mathworks.com/support

and search for "higher order ODE". The first two hits each include examples
of how to do that.

Alternately, you could use ODE15I, which can handle ODEs of the form f(t, y,
y') = 0.

--
Steve Lord
slord@mathworks.com

Thank you Torsten and Steve, Torsten your example works well, I just have to incorparate the inatial condtions, which I know how to do.

Thanks again.

Amir