Event function for ODE

5 次查看（过去 30 天）
José Anchieta de Jesus Filho
回答： Walter Roberson 2021-5-18
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

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2021-5-18
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.

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by