# Strange behavior with ODE45 dependent on tspan

조회 수: 2 (최근 30일)
Lorenzo Carletti 2024년 4월 12일
댓글: Bruno Luong 2024년 4월 13일
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이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Bruno Luong 2024년 4월 12일
편집: Bruno Luong 2024년 4월 12일
• 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

댓글을 달려면 로그인하십시오.

### 채택된 답변

Bruno Luong 2024년 4월 12일
편집: Bruno Luong 2024년 4월 12일
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이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Bruno Luong 2024년 4월 13일
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

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

R2023a

### Community Treasure Hunt

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

Start Hunting!

Translated by