Event function for ODE
5 views (last 30 days)
Show older comments
Hello guys, I would like your help with my code. I am trying to find the steady state of a mass-spring-damper system that is being excited by a sinusoidal force.
For this I imposed the condition that, when the difference between the peaks is close to 0, the integration should stop, using "Events", but I am not having success. I would greatly appreciate someone's help.
function EventFunctTest
clear
close all
clc
w0 = 2; % (Hz)
% IC
x0 = 0; v0 = 0;
IC = [x0,v0];
opts = odeset('Events',@events,'RelTol',1e-6,'AbsTol',1e-6);
[time,state_values2] = ode45(@(t,s) gg(t,s,w0),[0 inf],IC,opts);
theta1 = state_values2(:,1);
plot(time,theta1)
end
function sdot = gg(t,s,ww)
Meq = 0.086;
Ceq = 0.1664;
Keq = 166.3629;
Feq = 0.6385;
f1 = @(t,ww) Feq*sin(ww*t*2*pi); % Force
sdot = [s(2);
(f1(t,ww) - Ceq.*s(2) - Keq*s(1))/Meq];
end
function [check, ist, dire] = events(t,s)
y1 = findpeaks(s(1));
dire = [];
ist = 1;
check = abs(y1(end) - y1(end-1)) - 1e-8;
end
0 Comments
Answers (1)
Walter Roberson
on 18 May 2021
y1 = findpeaks(s(1));
s(1) is a scalar. findpeaks() is not going to do anything useful there.
The event function will not receive any information about the history of the ode. You could hypothetically record the history as it runs, but that is risky, as the ode will get evaluated at a number of points that will never be passed through in practice, as Runge-Kutta methods require sampling at nearby locations in order to estimate the local slopes.
If you need a reliable history of the function, you need to use one of the delay differential routines.
You do not have a constant delay (or at least you do not know what it is) so if you have any hope at all it would be in https://www.mathworks.com/help/matlab/math/state-dependent-delay-problem.html . This is not a topic I have enough experience with to know whether that can be usefully adapted.
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!