second order system of differential equations in Matlab

amodiitm

Hi every one !

I am a final year Engineering student of IIT Madras, India. I am doing a project(finite element analysis of a structure) which requires the solution of a system of second order differential equations. equation looks like below:

[M][U"]+K[U]=[F(t)]

M : Mass Matrix of size 202x202
K: Stiffness matrix of size 202x202
F(t): time dependent force matrix of size 202x1
M , K and F(t) matrices are known

U: displacement matrix for nodes , size 202x1
U'': double derivative of displacement matrix, size 202x1

Initial conditions :
[U]=0 at t=0
[U'']=0 at t=0

to find :

[U] at t

this is the displacement matrix for the nodes of the system of elements


I will be extremely happy if someone could come out with a MATLAB code to solve this problem, this will be a great help for me as I am struggling with this from past few days and not able to make any progress in the project. Thanks a lot.. eagerly waiting for reply..


Amod

D H

You can transform the second-order ODE to a first-order ODE by adding the velocities [u'] to your state. For lack of a better term I will denote this augmented state as x:

[tex]
\begin{aligned}
\mathbf x &= \bmatrix \mathbf u \\ \mathbf u' \endbmatrix \\
\mathbf x' &= \bmatrix \mathbf u' \\ \mathbf u'' \endbmatrix
\end{aligned}
[/tex]

In other words, your 202 element second order ODE becomes a 404 element first order ODE. Matlab is quite capable of handling such things.

amodiitm

hello,
Thanks for your reply,but I tried the same yesterday which did not yield any result, I will paste here the code

Fo=@ (t,i) F(i).*sin((2*pi()/T)*t);

Fo(1,2)
tspan=[0:100];
u0=zeros(202,1);
dudt0=zeros(202,1);
y0=[u0;dudt0];

[t,y]=ode45(@f,tspan,y0);

could you please give me the code if you can..

D H

I can't do that. That is against the rules of this forum. We have these rules because our goal at Physics Forums is to help you learn. You won't learn nearly as much if someone gives you the answer as compared to working it out yourself.

That said, why would you expect to get anything of interest from a system that starts in the equilibrium position and has an initial velocity of zero?

amodiitm

though the system may have initial velocity of zero, but there will be displacement of the nodes due to the exciting force that is on the right side of the equation which is a function of time which will cause displacement of the nodes with every passing time. I am interested in finding out the same displacement matrix after time t(value of this t is of our choice).

Actually, I have been trying this out from last few weeks and not able to solve, now i have very limited time to complete my project which will not move ahead without solving this therefore I asked for the code if possible..

D H

You said earlier "but I tried the same yesterday which did not yield any result". Surely you got some result. A compiler error is a kind of result. Not the result you want, but a result. A zero vector is also a result. What do you mean by "did not yield any result"?

Matlab has a good debugger. Have you tried using it?

amodiitm

sorry for sounding vague , actually excerpts from my code that i pasted in previous post, I suppose is based on the same principle that u suggested but i could not go ahead with it,as in there was problem in implementing it. I could not completely write the code itself..therefore i asked for your code if u could give.. nevertheless, I will try again the same and I should be able to tell you what problem I am facing exactly.. thanks a lot

amodiitm

Hi

I am little confused about how to implement the

x=[u;u'];

x'=[u';u''];

thing in actual code.

If I can solve two 2nd order differential equation in matrix form then I think I can write for the 202x202 matrix too.

But I could not write the exact code........ ( Hence I am not expecting any result..... matlab debugger is not of much help too)

Can you help me to write a code for just two 2nd order equation...

i.e

[1 2;3 4] [u'']+[4 6; 6 7] [u]=[4*sin(2*t) ;6*sin(2*t)]

Thanks again for your help.
Amod

D H

Matlab has a fairly nice debugger. I suggest you try it. Try putting a breakpoint in your derivative function. (You do have a derivative function, don't you?)

Rather than attacking the big problem, try a very simple one, such as a point particle undergoing constant acceleration. This is a very simple second-order ODE. Your state is a two vector, position and velocity. Call this state 'x': x(1)=position, x(2)=velocity. The state derivative is then x'(1)=x(2), x'(2)=a. Try writing a derivative function for this simple problem and use some integrator such as ode45 to propagate.

Karto

Hello!

First write in MATLAB

help assempde

I think this can solve your problem. If not, then write

help pde

This will give you a list of very usefull functions to solve PDE with MATLAB (but in 2-D only).

For a first approximation you can use a GUI called PDETOOL (so write pdetool in MATLAB).

After playing a bit with it, you can use assempde for stationary problems, parabolic for problems regarding the first derivative of the time, and so on.

I hope this will help you.

Kind regards.

John Creighto

The original poster converted it to an ODE. Thus perhaps help ODE might get him/her the right information. I’d consider using an ODE solver that works for stiff systems. In such a case the numeric solver should be given the Jacobin to speed up integration and improve accuracy.

trambolin

Hey, if you are still stuck at the point, try this!

Let [tex] u(t) = x_1(t), \dot u(t) = x_2(t) [/tex], and also implied by the previous ones [tex] x_2(t) = \dot x_1(t)[/tex]

[tex]x = \left[ \begin{array}{cc} x_1 &x_2 \end{array}\right]^T[/tex]

then the following matrix includes all the dynamics right? (Check this!)

[tex]\dot x = \left[ \begin{array}{cc} 0 &I \\ -(M^{-1}K) &0\end{array}\right]x+ \left[ \begin{array}{c} 0 &M^{-1}F(t) \end{array}\right][/tex]

Then, either use the initial command with your initial condition vector or in Simulink, create a State space form by choosing C as identity or which node displacement that you want to observe.