161 views (last 30 days)

Show older comments

Greetings !

I have obtained displacement time history of a sytem using ODE15i. I need to obtain FFT curve of this response. As ODE solvers use adaptive algorithm, the time interval between steps is not constant. I have tried both hanning windowing option and zero padding, but not successful. It is giving broad spectrum instead of peaks at certain frequencies. I tried NUFFT command too, but it doesn't seem to help it out. If I fix the time step while solving ODE, then only I get good FFT curve. I guess, It is not good option to fix the time step in adaptive algorithms. Here, I am attaching the code and the data for your reference which contain both fixed step and variable step data.

Kindly suggest what should be done in this scenario. Thank you.

Paul
on 24 Jul 2021

You can try using the tspan input to od15i to specify equally spaced sample points in time at which the solution is desired. Check

doc ode15i

for more information.

Paul
on 3 Sep 2021 at 0:30

1. I only used ode45 to illustrate that an explicit solver is feasible. Another solver (along with appropriate options) might be more appropriate for your problem. This link (Choose an ODE Solver) may be of interest if you haven't read it already.

2. Without seeing an actual example, it's hard to say anything more about resetting the initial conditions at the transitions, other than to say it sounds like you're holding y(1) and y(3) fixed, and allowing y(2), y(4) and yp(1:4) to float.

3. Here is my attempt to have ode15i return the second derivatives.

global t0 A1 A2 a B1 B2 C1 C2 D1 D2

% define parameters

A1 = 300;

A2 = 300;

a = 0.3*A1;

B1 = 1E6;

B2 = 1E6;

C1 = 1E2;

C2 = 1E2;

D1 = 1E3;

D2 = 1E3;

t0 = 0.1;

% initial condition for BASIC EQUATIONS

% y(1) and y(3) must be zero

y0=[0;0;0;0];

yp0=[0;0;0;0];

tspan_a = [0 1];

Even though the stated desire is to have y(1) and y(3) initialized to zero, the call to decic allows them to vary. The result of the call to decic is that y(1) ~=0. Maybe use the fixed_y0 input to decic to hold y(1) and y(3) at zero.

[y0_new,yp0_new,resrnm1] = decic(@diff,0,y0,[],yp0,[])

The call to ode15i should use y0_new as the initial condition as that was the output from decic().

%[t,y]=ode15i(@diff,tspan_a,y0,yp0_new,options); % this should be y0_new!

odefunz = @(t,z,zp) ([diff(t,z(1:4),zp(1:4)) ; z(5:6)-zp([2 4])]);

options = odeset('RelTol',1e-6,'AbsTol',1e-8,"Jacobian",@Fjac2,"Events",@dydt1Events);

Might want to consider setting MaxSep and InitialStep.

z0 = [y0_new; yp0_new(2); yp0_new(4)];

zp0 = [yp0_new; 0; 0];

[t,z] =ode15i(odefunz,tspan_a,z0,zp0,options);

plot(t,z(:,1:4)),grid

This result looks identical(ish) to the result in the comment above using ode45

Now check that z(5:6) returned from ode15i are the second derivatives of y(1) and y(3), or the first derivatives of y(2) and y(4), i..e,. yp(2) and yp(4)

subplot(211);plot(t,z(:,5),t,gradient(z(:,2),t),'o'),grid

subplot(212);plot(t,z(:,6),t,gradient(z(:,4),t),'o'),grid

Looks very reasonable.

function dydt = diff(t,y,yp)

global t0 A1 A2 a B1 B2 C1 C2 D1 D2

% define load

if t<=t0

F = 1E6*(1-t/t0);

else

F=0;

end

% second order coupled differential equations (BASIC EQUATION)

% yp(1) is 1st derivative of y1 and written as y(2) ;

% yp(2) is 2nd derivative of y1 ;

% yp(3) is 1st derivative of y3 and written as y(4) ;

% yp(4) is 2nd derivative of y3 ;

%% A1*d2/dt2(y1)+a*(d2/dt2(y1)-d2/dt2(y3))+B1*y1+D1*(y1-y3)+C1*(d/dt(y1)-d/dt(y3))=F

%% A2*d2/dt2(y3)+a*(d2/dt2(y3)-d2/dt2(y1))+B2*y3+D2*(y3-y1)+C1*(d/dt(y3)-d/dt(y1))=0

dydt=[yp(1)-y(2);

yp(3)-y(4);

A1*yp(2)+a*(yp(2)-yp(4))+B1*y(1)+D1*(y(1)-y(3))+C1*(y(2)-y(4))-F;

A2*yp(4)+a*(yp(4)-yp(2))+B2*y(3)+D2*(y(3)-y(1))+C2*(y(4)-y(2))];

end

% jacobian for (BASIC EQUATION)

function [dfdy, dfdp] = Fjac1(t,y,yp)

global A1 A2 a B1 B2 C1 C2 D1 D2

dfdy=[0 -1 0 0;

0 0 0 -1;

B1+D1 C1 -D1 -C1;

-D2 -C2 B2+D2 C2];

dfdp=[1 0 0 0;

0 0 1 0;

0 A1+a 0 -a;

0 -a 0 A2+a];

end

function [dfdz, dfdzp] = Fjac2(t,z,zp)

[dfdy, dfdyp] = Fjac1(t,z(1:4),zp(1:4));

dfdz = blkdiag(dfdy,eye(2));

dfdzp = [dfdyp zeros(4,2); -[0 1 0 0 0 0; 0 0 0 1 0 0] ];

end

function [value,isterminal,direction] = dydt1Events(t,y,yp)

value = y(2);

isterminal = 1;

direction = -1;

end

Mateo González Hurtado
on 25 Aug 2021 at 9:13

Hi,

If it is not your intention to recalculate the response with ODE: have you tried to interpolate the non-uniform data?

Regards,

S Priyadharshini
on 2 Sep 2021 at 9:24

doc ode15i

Walter Roberson
on 2 Sep 2021 at 9:41

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!