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

1 view (last 30 days)
Newbie
Newbie on 21 Sep 2021
Edited: Paul on 21 Sep 2021
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

Answers (1)

Paul
Paul on 21 Sep 2021
Edited: Paul on 21 Sep 2021
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

Community Treasure Hunt

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

Start Hunting!