Clear Filters
Clear Filters

Strange behavior with ODE45 dependent on tspan

6 views (last 30 days)
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
Bruno Luong on 12 Apr 2024
Edited: Bruno Luong on 12 Apr 2024
  • 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

Sign in to comment.

Accepted Answer

Bruno Luong
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
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

Sign in to comment.

More Answers (0)

Tags

Products


Release

R2023a

Community Treasure Hunt

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

Start Hunting!