Event function for ODE

5 ビュー (過去 30 日間)
José Anchieta de Jesus Filho
José Anchieta de Jesus Filho 2021 年 5 月 17 日
回答済み: 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 ExchangeOrdinary Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by