How to code jump of the solution of diff. eq. at discrete points?

I made a matlab program for the solution of impulsive differential equation. The solutions have jumps at the following sequence.

where \theta, N_0, \Gamma_a are contant values and 'l' is natural number. I need to give jump to the solution at the points of the above sequence. I am confused how to do. Please help me.

function impulse_diff
clc;
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 i = 1:n
 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
     %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]);
end

Answer 1

[Abhimenyu]

Hi Rakesh,

I understand that you want to code jumps of the solution of the impulsive differential equation. Discontinuities produce discrete approximations in the system, which can lead to inaccurate results since the events are not being counted after a jump. The jump has been coded correctly in the comments of the code provided.

%   if (mod(i,N0)==0) 
%    h = b; 
%     t(i+1) = t(i)+h; 
%   else 
%    h = theta; 
%     t(i+1) = t(i)+h; 
%   end 

This code can help to add the jump in the solution but the events after will have ‘Nan’ values because the solution will be stopped after the jump but the ‘for’ loop will continue for a fixed iteration number.

To effectively code a jump for solving an impulsive differential equation, ‘ode45’ function of MATLAB can be used. It is a custom implementation of the ‘Runge-Kutta’ method. The ‘ode45’ function is more efficient and accurate, and it also adapts the step size automatically according to the local error. Also, a ‘while’ loop will allow you to iterate over time steps flexibly and stop the loop when the final time is reached, instead of using a fixed number of iterations. 

Please refer the MATLAB code below that involves the discontinuity mentioned in the problem description and solves the differential equation to achieve results.

% function definition has been removed to visualise results. 
% The definition can be added back to make this a function. 

% step sizes 
t0 = 0; 
tfinal = 5; 

% Initialize jump parameters 
theta = 0.01; 
N0 = 2; 
Gamma_a = 0.7; 
b = N0 * (Gamma_a - theta) + theta; 

% initial conditions 
x = [-5; 3]; % vector of x1 and x2 
t = t0; 

% functions handles (used as arrays to reduce code length and increase readability) 
Fx = @(t, x) [-1/2 * x(1) - x(1).^3 + x(1) .* x(2).^2 - ((x(1).^2 + x(2).^2).^1/4) .* sign(x(1)); ... 
              -1/2 * x(2) - x(2).^3 - x(2) .* x(1).^2 - ((x(1).^2 + x(2).^2).^1/4) .* sign(x(2))]; % vector of Fx1 and Fx2 

% initialize arrays to store the values of t and x 
T = []; 
X = []; 

% loop until the final time is reached 
while t < tfinal 
    
    % append the current values of t and x to the arrays 
    T = [T; t]; 
    X = [X; x'];

    % update the step size according to the jump sequence 
    if (mod(length(T), N0) == 0) 
        h = b; 
    else 
        h = theta; 
    end 
    
    % solve the differential equation using ode45 with the current initial conditions and step size 
    [t_sol, x_sol] = ode45(Fx, [t, t + h], x);

    % update the values of t and x with the last values of the solution 
    t = t_sol(end); 
    x = x_sol(end, :)'; 
end 

% plot 
plot(T, X(:, 1)); 
hold on; 
plot(T, X(:, 2)); 
xlim([0, tfinal]); 
ylim([-6, 6]); 

This will help you to code the jumps and achieve a solution for the differential equation.

Please refer to the below MATLAB documentation link to know more about ‘ode45’ function:

https://www.mathworks.com/help/matlab/ref/ode45.html

I hope this helps to resolve your issue.

Thank you,

Abhimenyu.