Event function for ODE

4 visualizaciones (últimos 30 días)
José Anchieta de Jesus Filho
José Anchieta de Jesus Filho el 17 de Mayo de 2021
Respondida: Walter Roberson el 18 de Mayo de 2021
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

Respuestas (1)

Walter Roberson
Walter Roberson el 18 de Mayo de 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.

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by