Problem in Octave code of impulsive differential equations?

178 Views Asked by At

In the image, I have written about the impulsive differential equation. I want to code the solution in Octave with Runge-Kutta 4th order method:

I have tried it. But I am not able to code the discrete part of the system.

clc;
clear all;
close all;

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

%initial conditions
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 i = 1:n
  
  %updates time 
  t(i+1) = t(i)+h;
  
 
 %updates x1 and x2

 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);
 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
%plot
plot(t,x1)
hold on;
plot(t,x2)
xlim([0, tfinal]);
ylim([-6, 6]);  
1

There are 1 best solutions below

0
On

Since this will produce a discrete approximation to your system, you may need to use event detection for detecting the discrete events in your system, if not, possibly none of the events will be counted. This is a hybrid system.

You could define the events, and tell Matlab what they are. (https://se.mathworks.com/help/matlab/math/ode-event-location.html)

Alternatively: A much less accurate way, but probably simpler way, is to approximate this, is to check if $t_n < some \; t_k < t_n + \Delta t$, i.e. checking if you are passing one of these discrete points in every iteeration in the numerical method. So, define your dynamics as a proper 'function', and include this is an if-else statement inside.

EDIT: Perhaps better than the approach above: integrate with RK between the discrete events, apply the discrete event to the right time point, then continue with the continuous simulation. This is probably OK when you have deterministic discrete events happening at known times, and if it does not happen too often.