- You should devide omega by 2*pi to get the resonance period frequency (noy angular frequency) the( compute fs from that?
- You should sample correctly the pushort pulse.
- Your pulse is not Gaussion but the derivative of the Gaussian.. The total net force is 0 so it will resonate differently depending on the phase of the excitation

# Strange behavior with ODE45 dependent on tspan

6 views (last 30 days)

Show older comments

Lorenzo Carletti
on 12 Apr 2024

Commented: Bruno Luong
on 13 Apr 2024

I'm writing code to solve a simple spring mass damper system excited by a very short gaussian pulse using ode45. I've noticed that the solution varies wildly with what is given for the tspan argument. The code I've pasted uses a sample frequency of 3 times the natural frequency to generate tspan, and the resulting solution looks as expected. However if I change the sampling frequency (say to 4*omega_n), or the number of sample points N, the solution changes drastically to something widlly incorrect. I suspect the problem is in the tight timing of the gaussian pulse, as when I increase the width the solution remains stable across almost any tspan.

I'd like to know why this behavior occurs, why the one value of sf makes the solution work when no other value does, and how to prevent this problem in the future. The only reason I was able to get a valid solution was by guessing a checking different values of sf until something worked.

clear

close all

clc

% System parameters

m = 2; %(kg)

k = 1250; %(N/m)

c = 0.1; %(Ns/m)

gaussian_amp = 3; %(N)

width = 0.01; %(s)

omega_n = sqrt(k/m);

% sampling parameters

sf = 3*omega_n;

N = 5*sf;

tspan=[0:N-1]/sf;

parameters = {m,c,k,gaussian_amp,width};

[t,z] = ode45(@sDoF_solver,tspan,[0,0],[],parameters);

x = z(:,1);

v = z(:,2);

g = make_ddt_gaussian(gaussian_amp,width,1,t);

plot(t,x)

hold on

yyaxis right

plot(t,g)

function x_dot = sDoF_solver(t,x,para)

M = para{1};

C = para{2};

K = para{3};

force_amp = para{4};

width = para{5}

g = make_ddt_gaussian(force_amp,width,1,t); % using an arbitrary 1 second delay for nice looking plots

% x(1) = x

% x(2) = v

x_dot = [x(2); -C/M*x(2) - K/M*x(1) + g/M]

end

function g = make_ddt_gaussian(A,width,delay,t)

arg = -pi*((t-delay)/width).^2;

gau = A*exp(arg);

g = 4.1327*(t-delay)/width.*gau;

end

##### 1 Comment

Bruno Luong
on 12 Apr 2024

Edited: Bruno Luong
on 12 Apr 2024

### Accepted Answer

Bruno Luong
on 12 Apr 2024

Edited: Bruno Luong
on 12 Apr 2024

I would divide the time interval in three parts, as the typical time scales are different during the pulse and outside the pulse. ode45 seems to do strange stuffs with non uniform tspan.

There is a bunch of sentenses on the doc page that warns about the dependency of oder45 solution to tspan

% Scystem parameters

m = 2; %(kg)

k = 1250; %(N/m)

c = 0.1; %(Ns/m)

gaussian_amp = 3; %(N)

width = 0.01; %(s)

omega_n = sqrt(k/m);

% sampling parameters

fres = 1/(2*pi)*omega_n;

sf = 3*fres; % minimum Nyquist with 2*fres

Tend = 6; % the end of time interval of the ode intervals

pulsedelay = 1;

% We break the integration interval in three parts

halfpulsewindow = 3*width;

T = [0, pulsedelay-halfpulsewindow, pulsedelay+halfpulsewindow, Tend];

Ts = 1/sf; % desired sample period

dt = [Ts, width/10, Ts]; % Desire time step on each interval

parameters = {m,c,k,gaussian_amp,width,pulsedelay};

% Iterative solve in each interval

z0 = [0 0]; % Initial state

t = 0; % initial time

z = z0(:).';

for k=1:length(T)-1

Tspan_k = T([k k+1]);

N = ceil(diff(Tspan_k)/dt(k)) + 1;

Tspan_k = linspace(Tspan_k(1), Tspan_k(2), N);

[tk,zk] = ode45(@sDoF_solver,Tspan_k,z0,[],parameters);

z0 = zk(end,:); % Update the initial Cauchy condition for next iteration

% concateneate results

t = [t; tk(2:end)];

z = [z; zk(2:end,:)];

end

x = z(:,1);

v = z(:,2);

g = make_ddt_gaussian(gaussian_amp,width,pulsedelay,t);

plot(t,x,"-")

hold on

yyaxis right

plot(t,g)

legend('x position', 'Pulse')

function x_dot = sDoF_solver(t,x,para)

M = para{1};

C = para{2};

K = para{3};

force_amp = para{4};

width = para{5};

delay = para{6};

g = make_ddt_gaussian(force_amp,width,delay,t);

% x(1) = x

% x(2) = v

x_dot = [x(2); -C/M*x(2) - K/M*x(1) + g/M];

end

function g = make_ddt_gaussian(A,width,delay,t)

arg = -pi*((t-delay)/width).^2;

gau = A*exp(arg);

g = 4.1327*(t-delay)/width.*gau;

end

##### 1 Comment

Bruno Luong
on 13 Apr 2024

This code runs different sf and show now consistency of the results

% System parameters

m = 2; %(kg)

k = 1250; %(N/m)

c = 0.1; %(Ns/m)

gaussian_amp = 3; %(N)

width = 0.01; %(s)

omega_n = sqrt(k/m);

Tend = 6; % the end of time interval of the ode intervals

pulsedelay = 1;

% We break the integration interval in three parts

fres = omega_n/(2*pi);

halfpulsewindow = 3*width;

parameters = {m,c,k,gaussian_amp,width,pulsedelay};

upsamplingtab = [2 4 8 16 32 64]; % upsampling factors from resonance frequency

marker = 'o+*^v!<>xsdpe';

color = 'rgbcmyk';

h = zeros(1,length(upsamplingtab));

for j = 1:length(upsamplingtab)

upsampling = upsamplingtab(j);

% sampling parameters

sf = upsampling*fres;

% Window break points

T = [0, pulsedelay+halfpulsewindow*[-1,1], Tend];

Ts = 1/sf; % desired sample period

dt = [Ts, width/10, Ts]; % Desired time step on each interval

% Iterative solve in each subinterval

t = 0; % initial time

z0 = [0 0]; % Initial state at t=0

z = z0(:).';

for i = 1:length(T)-1

Tspan_i = T([i i+1]);

N = round(diff(Tspan_i)/dt(i)) + 1;

Tspan_i = linspace(Tspan_i(1), Tspan_i(2), N);

[tk, zk] = ode45(@sDoF_solver, Tspan_i, z0, [], parameters);

z0 = zk(end,:); % Update the initial Cauchy condition for next iteration

% concateneate results

t = [t; tk(2:end)]; %#ok

z = [z; zk(2:end,:)]; %#ok

end

x = z(:,1);

if upsampling == max(upsamplingtab)

linestyle = '-k';

else

linestyle = [marker(j) color(j)];

end

hold on

h(j) = plot(t,x, linestyle);

end

g = make_ddt_gaussian(gaussian_amp,width,pulsedelay,t);

yyaxis right

plot(t,g,'r','linewidth',1)

legend(h, arrayfun(@(x) sprintf('upsamplig=%g', x), upsamplingtab, 'unif', 0));

xlim([0.9 1.5])

function x_dot = sDoF_solver(t,x,para)

[M,C,K,force_amp,width,delay] = deal(para{:});

g = make_ddt_gaussian(force_amp,width,delay,t);

x_dot = [x(2); -C/M*x(2) - K/M*x(1) + g/M];

end

function g = make_ddt_gaussian(A,width,delay,t)

arg = -pi*((t-delay)/width).^2;

gau = A*exp(arg);

g = 4.1327*(t-delay)/width.*gau;

end

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!