Can someone help me implement my Heaviside function into my Kuramoto model please?

2 次查看(过去 30 天)
Hi there,
I have a Heaviside function that is 1 only when 5.76<t_pert<5.96 and I want to implement that in my Kuramoto model. However, this does not work as it tell me as my Heaviside is symbolic variable but my dde23 inputs are numerical. Please can anyone help?
Here is my code so far:
t1=5.76;
t2=5.96;
t_pert=linspace(0,20*pi,1000);
% Define the Heaviside function
syms t_pert
H=1-(heaviside(t_pert-t1)*heaviside(t2-t_pert));
fplot(H,[5.6,6])
title('Heaviside function')
N=2;
omega=2*pi;
tau=11*pi/36;
h=0.1; % step size
iter=100; % number of iterations
t=1:h:h*iter; % time
K_coeff=[10 10];
K_guess=[0 K_coeff(1,1);K_coeff(1,2) 0]; % full coupling matrix
tspan=[0 20*pi]; % interval of integration
s=ones(2,1); % history of DDEs
sol=dde23(@(t,theta,Z)(my_kuramoto_pert(t,theta,Z,omega,H,K_coeff)),tau,s,tspan)
function dthetadt_pert=my_kuramoto_pert(t,theta,Z,omega,H,K_coeff)
thetalag1=Z(:,1);
dthetadt_pert=[omega+H+K_coeff(1,1)*sin(theta(2)-theta(1));
omega+K_coeff(1,2)*sin(thetalag1(1)-theta(2))];
end

回答(1 个)

Paul
Paul 2021-9-21
编辑:Paul 2021-9-21
The function H seems to do the opposite of what you want it do, i.e., it's zero over the desired range and one outside that range:
t1=5.76;
t2=5.96;
% Define the Heaviside function
syms t_pert
H=1-(heaviside(t_pert-t1)*heaviside(t2-t_pert));
fplot(H,[5.6,6],'-o')
Once you have H defined the way you want, you can do this (the code executes, don't know if the result is what you expect):
Hfunc = matlabFunction(H); % convert H to a function that can be evaluated with numeric input
You could also define Hfunc directly without going through the symbolic stuff:
Hfunc1 = @(t)(double(t<5.76 | t > 5.96));
fplot(Hfunc1,[5.6 6],'-o')
N=2;
omega=2*pi;
tau=11*pi/36;
h=0.1; % step size
iter=100; % number of iterations
t=1:h:h*iter; % time
K_coeff=[10 10];
K_guess=[0 K_coeff(1,1);K_coeff(1,2) 0]; % full coupling matrix
tspan=[0 20*pi]; % interval of integration
s=ones(2,1); % history of DDEs
sol=dde23(@(t,theta,Z)(my_kuramoto_pert(t,theta,Z,omega,Hfunc,K_coeff)),tau,s,tspan); % use Hfunc instead of H
function dthetadt_pert=my_kuramoto_pert(t,theta,Z,omega,H,K_coeff)
thetalag1=Z(:,1);
dthetadt_pert=[omega+H(t)+K_coeff(1,1)*sin(theta(2)-theta(1)); % t is the input argument to H, which is really Hfunc
omega+K_coeff(1,2)*sin(thetalag1(1)-theta(2))];
end

类别

Help CenterFile Exchange 中查找有关 Assumptions 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by