Strange behavior with ODE45 dependent on tspan

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

[Bruno Luong]

• 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

Answer 1

[Bruno Luong]

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

[Bruno Luong]

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