15,961,997 members
1.00/5 (1 vote)
See more:
I have impulsive differential equation. I want to code them in octave.

What I have tried:

```clc;
clear all;
close all;

% step sizes
t0 = 0;
tfinal = 5;
h = 0.01;
n = ceil((tfinal-t0)/h)+1;

%initial conditions
s(1) = 0.5;
theta = 0.01;
N0 = 2;
Ta = 0.7;
b=N0*(Ta-theta)+theta;

x1(1) = -5;
x2(1) = 3;
t(1) =0;

% functions handles
Fx1 = @(t,x1,x2) -1/2*x1-x1.^3+x1.*x2.^2-((x1.^2+x2.^2).^1/4).*sign(x1);
Fx2 = @(t,x1,x2) -1/2*x2-x2.^3-x2.*x1.^2-((x1.^2+x2.^2).^1/4).*sign(x2);
for N = 1:n

if (mod(N,N0)==0)
s(N+1) = s(N)+b;
else
s(N+1) = theta;
end

for i = 1:ceil(s(N+1))
t(i+1) = t(i)+h;

##  if (mod(i,N0)==0)
##    h = b;
##    t(i+1) = t(i)+h;
##  else
##    h = theta;
##    t(i+1) = t(i)+h;
##  end
k1x1 = Fx1(t(i),     x1(i),       x2(i)        );
k1x2 = Fx2(t(i),     x1(i),       x2(i)        );
k2x1 = Fx1(t(i)+h/2, x1(i)+h/2*k1x1, x2(i)+h/2*k1x2);
k2x2 = Fx2(t(i)+h/2, x1(i)+h/2*k1x1, x2(i)+h/2*k1x2);
k3x1 = Fx1(t(i)+h/2, x1(i)+h/2*k2x1, x2(i)+h/2*k2x2);
k3x2 = Fx2(t(i)+h/2, x1(i)+h/2*k2x1, x2(i)+h/2*k2x2);
k4x1 = Fx1(t(i)+h,   x1(i)+h*k3x1, x2(i)+h*k3x2);
k4x2 = Fx2(t(i)+h,   x1(i)+h*k3x1, x2(i)+h*k3x2);

if (i== ceil(s(N+1)))
x1(i+1) = 0.5*x1(i)+h/6*(k1x1 +   2*k2x1 + 2*k3x1 +k4x1);
x2(i+1) = 0.8*x2(i)+h/6*(k1x2 +   2*k2x2 + 2*k3x2 +k4x2);
else
x1(i+1) = x1(i)+h/6*(k1x1 +   2*k2x1 + 2*k3x1 +k4x1);
x2(i+1) = x2(i)+h/6*(k1x2 +   2*k2x2 + 2*k3x2 +k4x2);
end
end
N = N+1;
end

%plot
plot(t,x1);
hold on;
plot(t,x2);
##xlim([0, tfinal]);
##ylim([-6, 6]);
drawnow();```
Posted
Updated 15-Jun-22 9:54am