how to add a perturbation while solving ode using ode45? can i use the analytical approximation of Heaviside function for it if yes how?

function dydt = rsf_ode4eq(t,y,d,par,v_lp,t_step,law)

dydt = zeros(4,1);
aa = par(1); bb = par(2); dc = par(3); sigma = par(4); stiff = par(5);
rad = par(6);B_s=par(7);d=par(8);m=par(9);r=par(10);p=par(11);
omega = y(1)*y(2)/dc;
dydt(1) = 1.d0 - omega;
%tau=1./(1+exp(-2*m*(t-d)));


% taudot =0.01*B_s*(2*m*exp(-2*m*(t-d))./((1+exp(-2*m*(t-d))).^2));
% if t <= t_step
%     v_loc = v_lp;
% elseif t<=t_step+d
% v_loc=1.00001*v_lp;
% else
%     v_loc = 1.00001*v_lp;
% end
if t <= t_step
    v_loc = v_lp;
elseif t <= t_step+d
    v_loc=1.000001*v_lp;
    
else
    v_loc = 1.00000*v_lp;
end

% taudot_step = sirac(t,d,m,B_s); 
% taudot_step=10e-8*B_s*(r-exp(-(t-d)/p));
% taudot_step=10-8*B_s*dirac(t-d);
%  x=dirac(t-d);
tau_dot  = stiff*(v_loc- y(2))+B_s*10e-6*x;
if strcmp(law,'A')
    dydt(1) = 1.d0 - omega;
elseif strcmp(law,'S')
    dydt(1) = -omega*log(omega);
end
dydt(2) = ((tau_dot/sigma - bb*dydt(1)/y(1)))/(aa/y(2)+rad/sigma);
dydt(3) = tau_dot;

end

solve:
clear
clc
close all

a = 0.021;
b = 0.025;
Dc = 100; % in microns
sigma = 1e1; % in MPa
mu_0 = 0.6;
m=10-9;
r=1;
p=1e-3;
Gmod = 3e4; % in MPa
b_s=mu_0*sigma;
v_shear = 3e9; % in microns/sec
rad = 0.5*Gmod/v_shear; % Radiation damping term
v_dyn = (a*sigma)/rad;
Kc = sigma*(b-a)/Dc; % in MPa/microns
v_init = 1e-2; % in microns/s
theta_init = Dc/v_init; % in seconds
tspan = [1e-5 3.091924568069068e+07]; % in seconds
solopt = odeset('RelTol',1e-16);
t_step = 1e5; % in seconds
Im_stiff_osc = sigma*(sqrt(b)-sqrt(a))^2/Dc;
rat_oscillatory_sol = Im_stiff_osc/Kc;
stiff = 0.9*Im_stiff_osc; % Stiffness 
tau_init=stiff*v_init;
B_s=mu_0*sigma;
% xlsread('t_0.xlsx');
% r_numf=xlsread('t_0','C2:C52');

law = 'A';
d=[2.5e7, 2.55e7, 2.6e7, 2.8e7 ,2.9e7 ,3e7];
par = [a,b,Dc,sigma,stiff,rad,B_s,d,m,r,p];

for i = 1:2%length(d)
    dd=d(i);
    [t,y] = ode45(@(t,y)rsf_ode4eq(t,y,dd,par,v_init,t_step,law),...
        tspan,[theta_init;v_init;tau_init;0],solopt);
    
    theta=y(:,1);
    v=y(:,2);
    
    semilogy(t,v,'b','linewidth',0.75);
    xlabel('Time[s]')
    ylabel('Velocity[m/s]')
    hold on
    grid on
    title('velocity vs time for perturbed system')
    tau= sigma*(mu_0 + a*log(y(:,2)/v_init) + b*log(y(:,1)/theta_init));
    hold on
    %     v_new(:,i)=v;
    %     t_new(:,i)=t;
    %     [pks,locs] = findpeaks(y(:,2),t);
    %     loc_new(:,i)=locs;
    
    yyaxis right
    semilogy(t,tau,'g','linewidth',0.75)
    hold on
    ylabel('\tau')
    xlabel('time[s]')
    title('\tau vs time[s] for perturbed system')
    %    
end 

[Star Strider]

It is not obvious to me what you want to do.

Likely the best approach to introducing a step input to your system is to intergrate it up to the time you want to introduce the step, save the last row of the ‘y’ ouitput (and the last ‘t’ value), add the step value to the appropriate column of that vector to use as the initial conditions of the subsequent integration, and then start it from that time to the desired integration total time.

Even the most robust numeric ordinary differential equation integrators cannot handle step inputs during the integration, since the step inputs are not differentiable.

[Sam Chak]

A total of 98 lines in the code without explanatory comments.

Where exactly do you want to add the perturbation?

dydt(1) = 1.d0 - omega;
dydt(2) = ((tau_dot/sigma - bb*dydt(1)/y(1)))/(aa/y(2)+rad/sigma);
dydt(3) = tau_dot;

[SOUVIK DARIPA]

to taudot @Sam Chak

Answer 1

[Sam Chak]

Hi @SOUVIK DARIPA

If the perturbation looks like a step disturbance, then you can try using the Logictic function instead, which is a continuously differentiable. In physical systems, there will be a slight delay between the swiching states. Thus, I think it is reasonable to use that. Three parameters describe how the graph looks like:

x     = linspace(-1, 2, 30001);
a     = 3;          % the supremum of the function
slope = 1e3;        % steepness of the curve (adjust as you wish)
xsp   = 0.5;        % the point at which the switching occurs
approx_step = a./(1 + exp(- slope*(x - xsp)));    % a logistic function
plot(x, approx_step, 'linewidth', 1.5), grid on
xlabel('x'), ylim([-1 4])

Then, try adding perturbation like this:

dydt(1) = 1.d0 - omega;

dydt(2) = ((tau_dot/sigma - bb*dydt(1)/y(1)))/(aa/y(2)+rad/sigma);

dydt(3) = tau_dot + approx_step;