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

4 views (last 30 days)
SOUVIK DARIPA on 1 Mar 2023
Answered: Sam Chak on 3 Mar 2023
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);
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(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
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;
law = 'A';
d=[2.5e7, 2.55e7, 2.6e7, 2.8e7 ,2.9e7 ,3e7];
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
SOUVIK DARIPA on 1 Mar 2023
to taudot @Sam Chak

Sam Chak on 3 Mar 2023
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(3) = tau_dot + approx_step;

### Categories

Find more on Mathematics in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by