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

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

[Sargondjani]

https://nl.mathworks.com/help/symbolic/conversion.html

[Newbie]

Thank you for the link, I tried using this but it does not work for me.

Answer 1

[Paul]

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